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An  investigation  of  the  turbulent  mixing  in  the  initial  region 
of  a  half-jet  composed  of  dissimilar  gas  streams  has  been  made.   An 
iso-energetic,  non-reacting  and  isobaric  system  with  various  velocity 
and  mass  flux  ratios  was  studied  both  experimentally  and  theoretically. 

The  flow  problem  was  formulated  as  an  initial  value  problem 
using  the  turbulent  boundary  layer  equations  in  conjunction  with 
phenomenological  models  for  the  turbulent  eddy  viscosity.   These 
models  consisted  of  Prandtl's  mixing  length  hypothesis,  Ferri's 
differential  mass  flux  formulation,  Schetz's  extension  of  the  Clauser 
integral  model  and  Alpinieri's  momentum  and  mass  flux  model.   The 
validity  of  the  above  models  for  the  turbulent  momentum  transport 
mechanisms  in  the  initial  region  of  a  half-jet  was   investigated  on 
a  consistent  basis.   The  eddy  diffusivity  was  obtained  from  the  eddy 
viscosity  model  by  considering  the  turbulent  Schmidt  number  as  a 
parameter.   The  analytical  solution  was  obtained  using  an  all-explicit 
forward  marching  finite  difference  scheme  with  the  stability  of  this 
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scheme  being  ensured  by  satisfying  the  von  Neumann  stability 
criterion. 

In  the  experimental  phase  of  the  study,  a  two-dimensional  test 
section  was  designed  and  built  to  be  used  in  conjunction  with  the 
existing  blow-down  wind  tunnel  facilities.   Interchangeable  nozzle 
blocks  were  designed  with  the  two-dimensional  method  of  characteristics 
to  provide  the  supersonic  flow.   A  primary  stream  of  air  at  Mach 
numbers  of  1.3  and  2.0  was  allowed  to  mix  with  a  secondary  stream  of 
either  helium  or  argon  at  a  Mach  number  of  1.3.   Quantitative  data  in 
the  form  of  total  and  static  pressure  measurements,  total  temperature 
and  mass  fraction  measurements  were  collected.   Qualitative  data  in 
the  form  of  schlieren  photographs  were  also  obtained.   These  data 
were  used  to  correlate  and  complement  the  analytical  study. 

Results  of  experimental  static  pressure  measurements,  both 
experimental  and  theoretical  excess  velocity  and  mass  fraction  profiles 
are  presented  in  graphical  form  for  all  configurations  tested. 


CHAPTER  I 
INTRODUCTION 

The  occurrence  of  free  viscous  layers  and  their  consequent 
effect  on  the  performance  of  many  contemporary  devices  has  stimulated 
considerable  study  of  these  flow  processes.   Some  examples  are  slot 
cooling,  supersonic  diffusion  flames,  dumping  of  fuel  at  high  velocities 
and  thrust  augmentation  in  jet  and  rocket  engines.   It  is  deemed 
worthwhile  to  briefly  describe  how  the  free  viscous  layers,  resulting 
from  the  mixing  of  cof lowing  streams,  are  encountered  in  some  of  the 
physical  applications  cited  above. 

1)  Slot  (or  film)  cooling:   The  possible  structural  failure 
of  modern  aircraft  due  to  excessive  heating  has  been  a  major  problem. 
This  problem  arises  in  rotating  machinery  components  of  aircraft 
engines  as  well  as  in  all  types  of  re-entry  and  high-speed  vehicles. 
Some  examples  are  turbine  blades,  rocket  nozzles  and  leading  edges  of 
hypersonic  aircraft.   One  method  of  solution  to  this  problem  lies  in 
slot  cooling.   Slot  cooling  is  a  method  whereby  a  coolant  gas  (or 
liquid)  is  injected  into  the  boundary  layer  of  the  surface  to  be  cooled. 
The  purpose  is  to  create  a  cold  film  layer  of  gas  between  the  surface 
and  the  hot  mainstream.   If  the  film  is  maintained  over  a  large  portion 
of  the  surface  to  be  cooled,  it  acts  as  a  partial  heat  shield  and 
therefore  reduces  the  heating  process  sufficiently  so  as  to  ensure 

safe  operation  of  the  vehicle. 

2)  Diffusion  flames:   A  mixing  controlled  combustion  has 
several  features.   The  heat  release  is  distributed  over  a  finite  length 
in  contrast  to  a  premixed  configuration  wherein  the  heat  may  be  released 


abruptly  through  a  detonation  process.   The  inherent  distribution  of 
heat  release  in  a  mixing  controlled  system  provides  a  mechanism  for 
obtaining  a  controlled  pressure  variation  enhancing  the  possible  use 
of  a  fixed  geometry  system.   Furthermore,  the  mixing  controlled 
combustion  process  can  take  place  in  supersonic  flow  eliminating  flow 
losses  (i.e.,  total  pressure  loss  due  to  shock  waves)  and  critical 
design  problems  required  for  subsonic  burning. 

3)   Dumping  of  fuel:   Mixing  of  an  injected  gaseous  fuel  in 
combustible  proportions,  although  desirable  in  combustion  chambers, 
is  usually  undesirable  when  fuel  is  vented  from  a  flight  vehicle.   In 
the  case  of  the  multistage  chemical  rocket  using  cryogenic  propellant 
such  as  hydrogen,  large  quantities  of  waste  gaseous  fuel  must  be 
vented  overboard.   After  the  waste  gas  is  dumped  overboard,  it  may  be 
exposed  to  regions  of  high  temperature  such  as  the  surface  boundary 
layer  and  the  nozzle  base  area  of  the  operating  first-stage  engine.  A 
solution  to  this  problem  is  the  venting  of  the  gas  in  such  a  manner 
that  the  mixture  is  diluted  below  the  lower  limit  of  flammability  in 
a  reasonably  short  distance  from  the  point  of  injection. 

A)   Thrust  augmentation:   In  addition  to  combustion  and  heat 
transfer  problems,  secondary  injection  is  of  importance  in  ejector 
systems  fcr  jet  and  rocket  engines.   High  velocity  secondary  gases  are 
ducted  into  the  nozzles  to  increase  the  total  exit  momentum  flux  of 
the  flow,  thus  obtaining  higher  thrust  vectors. 

For  any  kind  of  analytical  study  to  be  conducted  on  the  above 
physical  problems  there  is  a  need  to  know  the  rate  of  growth  of  the 
shear  layer,  subsequently  referred  to  as  the  mixing  region,  the 


variables  affecting  this  growth  and  the  velocity,  concentration  and 
temperature  profiles  inside  the  mixing  region. 

The  mixing  process  may  take  place  under  a  variety  of  conditions 
which  determine  the  method  of  solution  to  the  problem.   Some  of  these 
conditions  are  summarized  below: 

i.   The  transport  mechanism  governing  the  mixing  process  may 
be  laminar  or  turbulent, 
ii.   The  mixing  process  may  be  steady  or  unsteady, 
iii.   Mixing  may  be  isobaric  or  take  place  in  the  presence  of 
pressure  gradients, 
iv.   The  compositions  of  the  mixing  streams  may  be  similar 
or  dissimilar, 
v.   The  flow  may  be  compressible  or  incompressible, 
vi.   The  geometry  of  the  system  may  be  two-dimensional  or 
axi-symmetric. 
vii.   Chemical  reactions  may  take  place  in  the  mixing  region 
or  the  flow  could  be  frozen, 
viii.   The  mixing  may  be  isoenergetic  or  non-isoenergetic. 
ix.   The  mixing  streams  could  be  contained  by  walls  or  the 
streams  may  be  semi-infinite. 
In  most  physical  problems  many  of  the  complicating  factors  are  present 
but  some,  have  to  be  neglected  or  simplified  in  order  to  obtain 
analytical  solutions. 

A  brief  investigation  of  the  mixing  problem  first  reveals 
that  the  mixing  processes  are  almost  unanimously  turbulent.   Second, 
one  finds  that  at  least  one  of  the  mixing  streams  will  be  compressible. 
Thus,  the  mixing  process  of  interest  must  account  for  compressibility 


of  the  streams.   A  third  observation  one  might  note  is  that  the 

two  streams  are  very  likely  to  be  at  different  thermal  levels,  requiring 

an  assessment  of  these  effects  as  well. 

For  all  cases  of  mixing  previous  investigators  [1]   have 
concluded  that  for  the  purpose  of  analysis  the  mixing  process  may  be 
divided  into  two  regions.   The  first  is  the  developing  region  where 
the  velocity  profiles  are  non-similar,  and  second  the  asymptotic  (or 
fully  developed)  region  where  the  velocity  profiles  are  self-preserving. 
Figure  1  shows  a  schematic  of  the  mixing  region  resulting  from  the 
"contact"  of  two  uniform  parallel  semi-infinite  streams  which  are 
initially  separated  by  a  thin  splitter  plate.   The  rate  of  growth  • 
of  the  mixing  region  is  determined  by  the  turbulent  transport  mechanisms 
of  mass,  momentum  and  energy. 

Approach 

There  are  essentially  three  methods  in  use  today  by  which  the 
mixing  of  parallel  cof lowing  streams  are  studied.   These  consist  of 
two  methods  that  have  been  in  use  for  some  time,  and  a  method  that  has 
become  available  within  the  last  several  years  as  a  result  of  the  rapid 
development  and  use  of  high-speed  computers.   The  three  methods  are: 
1)  the  simple  momentum  integral  method,  2)  solution  of  the  equations 
of  motion  where  certain  assumptions  are  made  which  render  these 
equations  tractable  to  existing  analytical  techniques,  and  3)  step-by-step 
numerical  solution  of  the  equations  on  high-speed  digital  computers. 
All  of  these  methods  depend  ultimately  on  experimental  turbulent  mixing 


Numbers  in  brackets  designate  references. 


data  since  an  essential  part  of  all  these  techniques  is  the  specifica- 
tion of  the  local  turbulent  transport  coefficients  (eddy  viscosity, 
diffusivity,  conductivity)  of  the  particular  turbulent  flow  in  question 
at,  at  least,  one  general  location  in  the  mixing  region. 

In  practical  applications  one  is  sometimes  interested  in  the 
initial  stage  of  mixing,  where  the  upstream  velocity  has  an  important 
effect.   The  classical  similarity  solutions  for  the  shear  layer  cannot 
account  for  the  effect  of  the  upstream  boundary  layer  on  the  mixing 
process.   The  profile  similarity  assumption  limits  validity  of  the 
solution  to  the  region  past  the  developing  region  [2,3].   Solutions 
for  the  nonsimilar  problems  may  be  derived  using  the  integral  technique; 
however,  this  application  is  limited  in  practice  to  reasonably  smooth 
profiles.   Profiles  that  cannot  be  approximated  by  analytical 
expressions  such  as  the  step  function,  or  profiles  that  exhibit  large 
boundary  layer  deficits  lie  outside  the  scope  of  the  integral  method. 
In  some  of  the  analyses,  transformations  are  made  to  obtain  closed- 
form  solutions  that  severly  limit  the  variation  in  flow  variables.   In 
reference  [4]  a  closed-form  solution  for  this  problem  is  obtained,  but 
the  linearization  of  transformed  equations  that  is  employed  again  limits 
the  variation  in  flow  variables. 

Therefore,  the  analytical  approach  taken  in  this  investigation 
was  to  develop  a  solution  that  would  be  numerical  in  nature  and  permit 
initial  variations  of  density,  velocity  and  temperature  profiles,  and 
would  be  valid  close  to  the  flow  inlet  as  well  as  far  downstream.   The 
flow  field  was  treated  by  employing  the  von  Mises  transformation  to 
the  boundary  layer  equations  and  utilizing  several  different  hypothesized 
models  for  the  eddy  viscosity.   The  resulting  conservation  equations 
were  solved  numerically  with  an  explicit  type  "marching"  technique. 


As  stated  earlier,  experimental  data  are  essential  to  evaluate 
the  validity  of  any  formulated  turbulent  transport  coefficient.   Such 
data  on  the  mixing  of  compressible  streams  of  different  compositions 
are  available  in  literature  only  in  scarce  quantity.   The  main  body 
of  available  material  deals  with  the  fully  developed  flow  region  of 
axi-syir.r.etric  jtts.   Data  on  the  two-dimensional  mixing  in  the  initial 
region  are  practically  nonexistent. 

Thus,  an  experimental  investigation  was  done  to  provide  data 
for  the  two-dimensional  mixing  of  dissimilar  streams.   The  data  were 
used  to  complement  and  verify  the  analytical  approach  taken  in  this 
study. 

The  gas  dynamics  facilities  in  the  mechanical  engineering 
department  were  adapted  to  utilize  different  species  of  gas  for  the 
nixing  analysis.   The  high  pressure  air  supply  system  was  utilized  in 
the  primary  stream  and  a  bank  of  commercial  bottled  gases,  i.e.,  argon 
and  helium,  was  used  to  supply  the  secondary  stream.   A  test  section 
was  designed  which  permitted  the  use  of  interchangeable  nozzle  blocks. 
The.  isoenergetic  mixing  took  place  in  this  two-dimensional  test  section 
where  both  the  primary  and  secondary  streams  were  supersonic.   Static 
and  total  pressure  measurements  were  made  in  the  test  section.   Gas 
samples  were  withdrawn  from  the  mixing  region  and  collected  in  a  series 
of  vacuum  bottles.   These  samples  were  subsequently  analyzed  with  a 
gas  chroma to graph.   A  schlieren  system  was  also  employed  to  obtain 
qualitative  results. 


CHAPTER  II 
LITERATURE  REVIEW 

The  analyses  of  the  mixing  of  turbulent  flow  fields  have  been 
performed  by  employing  the  laminar  flow  boundary  layer  equations  modified 
by  replacing  the  laminar  viscosity  with  the  eddy  viscosity  and  by 
replacing  the  laminar  Prandtl,  Lewis  and  Schmidt  numbers  by  their 
turbulent  counterparts.   The  use  of  the  boundary  layer  equations  is 
justified  by  the  fact  that  the  region  of  space  in  which  a  solution  is 
being  sought  does  not  extend  far  in  the  transverse  direction,  as 
compared  with  the  main  direction  of  flow,  and  that  the  transverse 
gradients  are  large.   The  assumptions  and  details  involved  in  the 
reduction  of  the  general  Navier-Stokes  equations  to  the  boundary  layer 
form  are  well  documented  in  references  [1],  1.2]  and  [3],  and  will  not 
be  reviewed  here.   The  two-dimensional  continuity  and  momentum  equations 
for  steady  flow  are  presented  below  to  aid  in  understanding  some  of 
the  assumptions  and  techniques  used  by  various  investigators  in 
obtaining  solutions.   These  equations  are  discussed  in  detail  in 
Chapter  III. 
Continuity:  3pu  +  3pv   n  (2  ]  ) 
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Y-Momentum:  C  «=  ~-  (2.3) 

(Here  x  denotes  the  turbulent  shearing  stress  cr  the  Reynolds  stress 
and  is  usually  expressed  in  terms  of  the  time  mean  average  of  the 
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The  works  of  previous  investigators  will  be  reviewed  in  three 
sections;  1)  nixing  of  semi-infinite  streams,  2)  mixing  of  contained 
(ducted)  streams,  and  3)  some  experimental  investigations  on  the 
mixing  of  dissimilar  gases.   The  first  two  sections  will  primarily  be 
involved  in  the  analytical  approaches.   The  last  section  is  included 
because  it  has  a  definite  bearing  on  the  experimental  aspect  of  the 
present  investigation. 

Mixing  of  Seni-Inf inite  Streams 

An  excellent  bibliography  of  both  experimental  and  theoretical 
work  on  turbulent  mixing  prior  to  1950  is  given  in  a  paper  by  Forstall 
and  Shapiro  [5].   Two  of  the  works  in  this  paper  should  be  mentioned 
here  since  they  formed  the  starting  point  and  basis  for  some  of  the 
more  recent  investigations. 

The  mixing  of  semi-infinite  incompressible  streams  was  considered 
analytically  as  early  as  1926  by  Tollmien  [6].   Using  a  similarity 
transformation  of  the  type  n  =  y/x  and  Prandtl's  [7]  mixing  length 
hypothesis   for  the  turbulent  transport  mechanism,  he  obtained  a 
numerical  solution  for  the  fully  developed  region  of  a  two-dimensional 
turbulent  jet  exhausting  into  a  quiescent  atmosphere.   This  solution 
was  later  extended  by  Kuethe   [8]  for  various  boundary  conditions. 

Gortler  [9]  utilized  Prandtl's  [10]  second  hypothesis  for  the 
eddy  viscosity  to  obtain  a  new  analytical  solution  for  the  incompressible 
mixing  of  two  parallel  streams.   The  result  was  a  series  solution  in 
contrast  to  Tollmien' s  numerical  solution  and  offered  the  further 


A  detailed  discussion  on  the  turbulent  transport  mechanisms 
is  presented  in  Chapter  IV. 


advantage  that  for  sufficiently  large  secondary  velocities,  the  velocity 
profile  could  be  approximated  by  the  error  function. 

Both  Tollmien's  and  Gortler's  solutions  involved  the  utilization 
of  a  stream  function  which  was  proportional  to  a  function,  F,  of  the 
similarity  variable  n  =  oy/x.   Thus,  the  partial  differential 
equations  of  motion  were  reduced  to  a  single  ordinary  differential 
equation: 

TcIJmien's  problem:      F"'  +  F  =  0 
Gortler's  problem:       F' ' '  +  2aFF''  =  0 

Tollmien  experimentally  determined  a  value  of  12  for  o   from  low  subsonic 
turbulent  mixing  experiments. 

It  was  desired  to  extend  this  type  of  solution  to  compressible 
flows.   Thus,  mixing  analyses  have  primarily  been  concerned  with  the 
development  of  the  theoretical  expressions  for  the  mixing  similarity 
parameters;  specifically,  the  evaluation  of  the  "spread  rate  parameter", 
a.   One  widely  used  method  is  to  apply  a  coordinate  transformation  to 
the  compressible  flow  equations  such  that  the  transformed  equations  are 
in  the  form  of  the  incompressible  ones.   Then,  data  from  both  the 
compressible  and  incompressible  domain  may  be  used  to  predict  a.      One 
such  method  is  the  extension  of  Howarth's  [11]  transformation  which 
was  developed  for  the  analysis  of  compressible  laminar  boundary  layers. 
The  transformation  is  essentially  stretching  of  the  y-coordinate  and  may 
be  defined  in  one  form  by 

* 

(2.4) 
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These  relations  will  reduce  the  compressible  flow  equations  to  the 
incompressible  flew  equations  in  the  absence  of  external  pressure 
gradients.  Mager  [12]  removed  the  restriction  on  the  pressure  gradient 
by  postulating  that  the  shear  stress  was  invariant  under  the  transfor- 
mation. 

There  are  a  number  of  theoretical  estimates  of  the  effects  of 

compressibility  on  the  spread  rate  parameter.   These  estimates  relate 

*  1  * 

the  ratio  a/c     to  the  free  stream  Mach  or  Crocco  numbers;  a  is  the 

similarity  parameter  for  incompressible,  isoenergetic  flow.   Abramovich 

[1]  estimated  lateral  turbulent  transport  using  the  classical  Prandtl 

hypothesis  for  eddy  viscosity  in  conjunction  with  a  hypothetical 

characteristic  longitudinal  velocity  in  the  shear  layer  in  order  to 

predict  the  growth  rate  of  the  mixing  region.   Bauer  [13]  based  his 

model  on  a  compressible  analog  of  the  mixing  length  concept;  together 

with  the  error  function  approximation  of  the  velocity  profile,  he 

was  able  to  estimate  the  spread  rate. 

Channapragada  and  Uoolley [14]  ,  using  the  Howarth  transformation 

in  conjunction  with  Mager 's  postulate,  reduced  the  governing  equations 

to  the  form  of  Tollmien's  or  Gortler's  problem  depending  on  the  model 

for  the  eddy  viscosity.   From  the  transformation  they  concluded  that 

the  parameter  a   varied  across  the  mixing  region  for  compressible  flow 

fields.   For  the  two  stream  mixing  problem,  they  related  a   to  the  total 

temperature  and  velocity  ratios  of  the  two  streams  in  addition  to  the 

primary  stream  Crocco  number.   This  model  for  o  agreed  well  with  the 


The  Mach  and  Crocco  numbers  are  related  through  the  equation 
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empirical  relation  of  Korst  and  Tripp  [15]  for  total  temperature  ratios 
of  unity,  and  exhibited  similar  trends  to  the  predictions  of 
Channapragada  [16].   This  work  was  later  extended  by  Wool  lay [17]  to 
the  case  of  two  dissimilar  streams  (with  same  specific  heat  ratios)  in 
the  presence  of  small  pressure  gradients. 

Laufer  [13]  applied  a  Howarth  type  of  transformation  to  the 
time-dependent  rather  than  the  mean  equations  of  motion,  as  had  been 
done  in  the  past.   When  the  transformed  equations  were  averaged  and 
the  correlation  between  the  fluctuations  of  the  temperature  and  velocity 
gradients  neglected,  the  incompressible  turbulent  equations  for  the 
mean  flow  were  obtained;  and  these  could  be  solved  with  conventional 
methods. 

Utilizing  Prandtl's  second  hypothesis  and  an  extension  of 
Warren's  [IS]  momentum  integral  method,  Donaldson  and  Gray  [20] 
analyzed  the  turbulent  mixing  and  decay  of  axially  symmetric,  compressible 
free  jets  of  dissimilar  gases.   They  concluded  from  a  comparison  of 
data  with  theoretical  results  that  a  general  relationship  existed,  at 
each  axial  position  in  the  jet,  between  a  local  mixing  rate  parameter 
and  the  local  Mac.h  number.   Furthermore,  this  relationship  was 
independent  of  the  physical  properties  or  the  thermodynamic  state  of 
the  mixing  gases  (i.e.,  independent  of  molecular  weight  and  enthalpy). 
This  method  of  analysis  was  later  extended  to  the  case  of  coflowing 
streams  by  Smoot  and  Purcell  [21]. 

Peters  [22]  presented  a  detailed  discussion  on  various  eddy 
viscosity  theories  and  compressibility  transformations  during  the 
development  of  a  transport  model  incorporating  a  dual  scale  of  eddy 
sizes  which  was  used  to  predict  o  in  the  compressible  regime.   Lamb  [23] 
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developed  a  theory  which  permitted  the  estimation  of  the  effect  of 
compressibility  and  heat  transfer  on  the  spread  rate  parameter  for 
fully  developed  mixing  zones.   Integral  forms  of  the  conservation 
equations  were  used  to  specify  the  flow  characteristics  along  the 
"dividing  streamline"  between  the  two  streams.   By  the  application  of 
the  Javier-Stokes  equations  to  this  streamline  he  was  able  to  calculate 
a  position  parameter  which  in  turn  yielded  an  expression  for  a. 

In  a  later  work  by  Lamb  and  Bass  [24]  an  analysis  of  the  methods 
involved  in  the  correlation  of  the  parameter  a  was  made.   It  was 
observed  that  differences  in  the  various  predictions  at  high  Mach 
numbers  (on  the  order  to  6-8)  were  as  much  as  50  percent.   It  was 
also  seen  that  the  trends  of  the  Channapragada  and  Woolley  theory 
appeared  to  be  opposite  to  those  of  the  other  analyses,  showing  a  large 
effect  of  compressibility  at  low  Mach  numbers  and  very  little 
influence  at  higher  values  of  M^. 

In  an  effort  to  compare  and  consolidate  different  theoretical 
velocity  profiles  and  expressions  for  a,  Korst  and  Chow  [25]  pointed 
out  that  the  spread  rate  parameter  depended  on   1)  the  selected  eddy 
viscosity  model,  2)  the  methods  of  theoretical  analysis  and  3)  the 
definition  of  profile  matching.   On  this  basis  they  established 
theoretical  relations  which  attempted  to  reconcile  discrepancies 
between  different  analytical  solutions  so  that  all  available  information 
on  a    could  be  utilized. 

It  can  be  concluded  that  many  of  the  compressible  mixing  investi- 
gations have  placed  an  emphasis  on  obtaining  similarity  solutions 
which  are  valid  for  regions  far  removed  from  the  initial  point  of 
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contact  of  the  two  streams.   Major  effort  has  gone  into  methods  which 
attempted  to  predict  the  behavior  of  the  spread  rate  parameter  in  the 
compressible  domain.   Unfortunately  there  is  no  widely  accepted  relation 
for  o.   This  may  be  coupled  to  the  fact  that  there  is  no  universal 
model  for  the  turbulent  transport  mechanism  due  to  the  lack  of 
understanding  of  the  physics  of  the  phenomena.   Another  drawback  is 
that  no  information  can  be  extracted  about  the  initial  and  transitional 
regions  of  mixing  from  the  spread  rate  parameter  or  similarity  solutions. 

Little  work  has  been  done  on  the  effect  of  nonuniform  initial 
velocity  profiles  on  the  mixing  process.   Wygnanski  and  Hawaleshka 
[26]  attempted  a  series  solution  to  the  case  of  the  incompressible 
asymmetric  jet  (i.e.,  created  by  the  mixing  of  a  wall  jet  with  quiescent 
surrounding  fluid  downstream  of  the  trailing  edge).   However,  the  region 
of  validity  of  this  method  as  well  as  its  accuracy  depended  critically 
on  an  accurate  knowledge  of  the  initial  velocity  profile  and  its 
derivatives  with  respect  to  y  and  the  number  of  terms  retained  in  the 
series  expansion.   Korst  and  Chow  [27]  used  integral  methods  to  account 
for  the  effects  of  initially  disturbed  profiles  on  the  mixing  region. 
Lamb's  [23,28]  "dividing  streamline"  solution  using  the  momentum 
integral  method  yielded  some  information  on  the  transition  region. 

Most  of  the  investigators  have  assumed  a  value  of  unity  for 
the  turbulent  Prandtl  and  Lewis  numbers.   This  assumption  enabled 
them  to  obtain  a  solution  to  the  momentum  equation  and  use  this 
solution  in  conjunction  with  the  Crocco  integral  relation  to  automatically 
satisfy  the  species  and  energy  equations.   Although  a  value  of  unity 
for  the  turbulent  Prandtl  number  may  be  justified  in  most  cases,  the 
turbulent  Schmidt  number  (Fr/Le)  has  been  shown  experimentally  to  vary 
between  0.5  and  2.0  [29,30,31]. 
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Another  widely  used  method  in  obtaining  solutions  to  the 
mixing  problem  is  one  which  is  based  on  the  linearization  of  the 
conservation  equations  in  the  plane  of  the  vcn  Mises  variables  x^hile 
retaining  the  essential  non-linear  nature  of  the  equations  in  the  physical 
plane.   This  results  in  a  linear  partial  differential  equation  in  the 
form  of  the  unsteady  heat  equation  which  can  be  treated  with  conventional 
methods  [32]  subject  to  initial  and  boundary  conditions.   However, 
relating  the  intrinsic  coordinate  system  to  the  physical  plane  by 
means  of  integral  forms  of  the  conservation  equations  is  rather  involved 
and  requires  the  specification  of  a  compressible  eddy  viscosity  model. 
The  works  of  Libby  [33],  Ferri  et  al.  [34],  Alpinieri  [35],  Kleinstein 
[4j,  and  Schetz  [36]  were  all  based  on  this  type  of  a  solution  with 
each  solution  deviating  from  the  other  in  the  formulation  of  the 
turbulent  transport  mechanisms.   These  formulations  will  be  analyzed 
in  Chapter  IV. 

Numerical  solutions  have  been  utilized  recently. to  obtain 
mixing  characteristics  throughout  the  flow  field  [37,38,39].   These, 
however,  required  appropriate  models  for  the  Reynolds  transport  terms 
for  each  region  of  mixing;  that  is,  the  initial,  transition  and 
fully  developed  regions. 

Mixing  of  Confined  (Ducted)  Streams 

Integral  techniques  have  been  used  in  general  to  obtain 
solutions  to  confined  turbulent  jet  mixing  problems.   In  most  of  the 
analyses  the  axial  velocity  and  the  shear  stress  were  assumed  to  obey 
similarity  laws,  so  that  expressions  for  the  turbulent  transport 
coefficients  were  not  needed.   In  1955,  Craya  and  Curtet  [40] 
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established  an  approximate  theory  for  confined  jet  mixing  of  streams 
of  identical  composition.   This  theory  was  further  developed  by  Curtet 
[41,42]  and  was  followed  by  additional  theoretical  and  experimental 
studies  by  Curtet  and  Ricou  [43]  and  Curtet  and  Barchilon  [44].   The 
theoretical  analysis  was  based  on  assumptions  of  zero  radial  pressure 
gradient,  uniform  and  non-turbulent  axial  velocity  outside  the  mixing 
region,  and  similarity  of  the  axial  velocity  profiles  inside  the  mixing 
region.   Experimental  observations  led  to  the  assumption  of  similar 
Gaussian  velocity  profiles  in  the  developing  region. 

In  1965,  Hill  [45,46]  carried  out  analytical  studies  of  an 
isothermal  homogeneous  confined  jet  mixing  in  order  to  predict  the 
mean  velocity  field  in  the  flow.   In  this  analysis,  an  integral 
technique  was  aJso  used  and  the  shear  integrals  were  evaluated  using 
free  jet  data.   However,  the  assumptions  made  were  such  that  the 
effects  of  a  confining  wall  were  not  significantly  taken  into 
consideration. 

Dealy  [47]  studied  the  effects  of  conditions  at  the  inlet  on 
the  flow  phenomena  in  a  confined  jet  mixing  system.   Dealy  concluded 
that,  for  systems  with  low  jet-to-confining  tube  radius  ratio,  the  flow 
in  the  near  regime  was  indeed  independent  of  the  nature  of  the  jet 
source;  had  similar  velocity  profiles  and  was  amenable  to  analysis 
by  the  common  momentum  integral  technique.   But  for  large  jet-to- 
confining  tube  radius  ratios,  the  mixing  mechanism  was  found  to  be 
dependent  strongly  on  the  flow  conditions  in  the  jet  exit.   Further 
experiments  by  Dealy  also  showed  that,  for  fully  developed  turbulent 
flow  at  the  jet  source,  mixing  took  place  more  rapidly  (because  of 
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larger  turbulent  stresses)  than  for  the  case  corresponding  to  a  uniform 
flow  at  the  jet  source. 

Trapani  [48]  carried  out  an  experimental  study  of  turbulent 
jets  with  solid  boundaries  in  the  transverse  direction  in  order  to 
investigate  their  application  in  certain  fluidic  devices.   Comparison 
with  the  flew  characteristics  of  a  two-dimensional  turbulent  free  jet 
showed  that  the  presence  of  solid  transverse  boundaries  definitely 
alters  the  behavior  of  the  flow.   The  bounded  jet  (i.e.,  the  jet 
bounded  by  plates  above  and  below)  was  seen  to  spread  less  rapidly 
than  the  free  jet.   On  the  other  hand,  the  confined  jet  (i.e.,  the  jet 
enclosed  on  all  sides)  was  observed  to  spread  more  rapidly  than  the 
free  jet;  this  effect  was  attributed  to  the  development  of  an  adverse 
axial  pressure  gradient  in  the  confined  flow. 

An  extensive  study  of  the  ducted  turbulent  mixing  process 
for  supersonic  flows  was  carried  out  by  Peters  et  al.  [39,49,50] 
experimentally  as  well  as  analytically.   An  integral  theory  for  the 
ducted  flow  was  presented  for  arbitrary  axi-symmetric  duct  geometry. 
Cases  of  both  frozen  and  equilibrium  chemistry  were  considered  as  the 
mode  of  chemical  reaction  in  the  mixing  zone.   At  initiation  of  mixing, 
the  boundary  layer  was  considered  negligible  as  were  the  viscous  effects 
at  the  duct  wall.   The  velocity  profiles  in  the  turbulent  mixing  zone 
were  assumed  to  be  similar  and  were  represented  by  a  cosine  function. 
The  turbulent  shear  stress  in  this  variable  density  mixing  layer  was 
treated  by  the  use  of  a  modified  Prandtl  eddy  viscosity  model.   The 
free  mixing  concept  of  shear  and  velocity  profile  similarity  were 
assumed  to  be  applicable  in  the  main  region.   The  turbulent  Prandtl 
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and  Lewis  numbers  of  unity  were  used  in  the  analysis.   From  the  results 
of  analysis  and  experiment,  it  was  concluded  that  the  integral  method 
developed  permitted  reasonably  accurate  computations  of  the  flow  in 
complex  mixing  systems  such  as  air-air  ejectors  and  air-augmented  rockets, 

Emmons  [51]  also  developed  an  analysis  for  predicting  the  flow 
characteristics  in  the  mixing  region  of  a  particle-laden  turbulent 
rocket  exhaust  and  the  surrounding  air  stream.   Neglecting  the  boundary 
layer  at  the  confining  wall,  the  turbulent  boundary  layer  equations 
were  used  to  describe  the  flow  in  the  mixing  zone.   The  eddy  viscosity 
model  was  assumed  to  vary  with  the  streamwise  coordinate.   The 
system  of  partial  differential  equations  governing  the  flow  was 
transformed  using  the  von  Mises  transformation  and  then  solved  by 
finite  difference  methods.   Similar  approaches  were  taken  by  Cohen 
[52],  Edelman  and  Fortune  [38]  and  Ghia   et  al.  [53]. 

Mixing  of  Dissimilar  Gases 

Major  Experimental  Investigations 

Perhaps  the  earliest  experimental  study  on  the  relative  rates 
of  diffusion  of  momentum  and  energy  in  turbulent  jets  was  done  by 
Ruden  [54],  who  found  that  energy  diffused  mere  rapidly  than  momentum 
in  isobaric  incompressible  jets.   Forstall  and  Shapiro  [5]  studied  the 
diffusion  of  mass  and  momentum  between  coaxial  jets  of  air  where  the 
central  jet  was  composed  of  approximately  10  percent  of  helium  for  use 
as  a  tracer  gas.   The  velocities  considered  in  their  experiments  were 
in  the  low  subsonic  range.   They  concluded  that  mass  diffusion  was 
more  rapid  than  momentum  diffusion.   Furthermore,  the  Schmidt  number, 
which  measures  the  relative  rates  of  transfer  of  mass  arid  momentum, 
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was  found  to  be  independent  of  the  velocity  ratio.   A  similar  result 
was  obtained  by  Keagy  and  Weller  [55]  from  their  experiments  on  helium 
and  carbon  dioxide  jets  exhausting  into  quiescent  air. 

Corrsin  and  Uberoi  [56],  using  a  heated  jet  exhausting  into 
a  quiescent  region  of  different  density,  investigated  the  effect  of 
large  density  differences  in  the  mixing  zone.   They  found  that  a  decrease 
of  jet  density  with  respect  to  that  of  the  receiving  medium  caused  an 
increase  in  the  rate  of  spread  of  the  jet. 

Isoenergetic  mixing  between  carbon  dioxide  and  hydrogen  central 
jets  exhausting  into  a  moving  concentric  stream  of  air  was  investigated 
by  Alpinieri  [35].   Using  velocities  in  the  low  to  high  subsonic  range, 
he  concluded  that  the  Schmidt  number  was  essentially  constant  and,  in 
agreement  with  previous  results,  that  mass  appeared  to  diffuse  more 
readily  than  momentum.   It  was  also  observed  that  no  tendency  toward 
segregation  of  the  two  jets  was  evident  when  either  the  velocity  ratio, 
mass  flow  ratio  or  momentum  flux  ratio  was  made  equal  to  unity.   Yates 
[57]  studied  the  supersonic  slot  injection  of  hydrogen  into  a 
supersonic  air  stream.   He  found  that  the  concentration  profiles  exhibited 
similarity  before  the  velocity  profiles  did.   The  rate  of  growth  of  the 
energy  and  concentration  layers  was  observed  to  be  about  the  same  and 
exceeded  the  growth  of  the  momentum  layer. 

Zakkay  e.t  al.  [58]  undertook  an  extensive  experimental  inves- 
tigation of  the  turbulent  mixing  of  two  dissimilar  gases.   The  axi- 
symmetric  mixing  analysis  was  carried  out  to  determine  the  turbulent 
transport  coefficients  for  hydrogen-,  helium-,  and  argon-air  mixtures. 
The  external  stream  of  air  was  maintained  at  a  constant  Mach  number  of 
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1.6  and  the  inner  jet  was  either  subsonic  or  supersonic.   Their 
conclusions  may  be  summarized  as,  1)  the  centerline  decay  was  not 
influenced  by  the  molecular  weight  or  the  initial  boundary  layer  of 
the  jet;  2)  the  radial  velocity  profiles  exhibited  similarity  past  the 
potential  core  and  could  be  correlated  by  the  mass  flux  ratio;  3)  no 
dependence  of  Schmidt  number  on  molecular  weight  could  be  observed; 
4)  the  deviation  of  the  turbulent  Schmidt  number  from  unity  was 
considerable  in  several  cases. 

However,  from  their  experimental  investigation  of  the  axi- 
symmetric  turbulent  jets  of  air,  helium  and  Freon  12,  Abramovich  et  al, 
[59]  concluded  that  it  was  ineffective  to  attempt  to  describe  the 
characteristics  of  the  jet  by  any  cojoint  complex  of  the  variables 
(p,u)  such  as  the  mass  flux  ratio  or  the  momentum  flux  ratio. 


CHAPTER  III 
ANALYSIS  FOR  THE  MIXING  OF  DISSIMILAR  STREAMS 

Introduction 


The  literature  survey  presented  shows  that  several  simplifying 
assumptions  were  made  to  render  the  analytical  model  tractable  from 
the  mathematical  point  of  view.   The  "boundary  layer"  forms  of  the 
conservation  equations  have  been  shown  to  be  applicable  to  the  cases 
of  both  confined  and  unconfined  mixing  problems.   In  a  great  majority 
of  the  investigations,  the  correlation  of  the  analysis  with  data  was 
restricted  to  the  main  region  of  mixing,  where  similarity  of  the 
velocity  profiles  was  assumed. 

In  the  present  study,  the  mixing  problem  is  formulated  as  an 
initial  value  problem  in  the  von  Mises  plane  using  the  turbulent 
boundary  layer  equations.   Unlike  the  laminar  problem,  the  transport 
mechanisms  in  turbulent  mixing  do  not  depend  on  the  fluid  properties 
alone,  but  also  on  geometric  and  dynamic  factors  of  the  flow  system. 
Determination  of  the  necessary  transfer  coefficients  by  an  exact 
theoretical  analysis  is  presently  not  possible  due  to  the  lack  of 
understanding  of  the  turbulence  phenomena;  therefore,  a  semi- 
empirical  approach  is  employed  in  order  to  correlate  the  data  obtained 
in  the  experimental  phase  of  the  investigation.   The  governing 
conservation  equations  of  the  flow  are  approximated  by  their  finite 
difference  forms  and  the  solution  is  obtained  by  an  explicit  numerical 
scheme.   Numerical  stability  is  ensured  by  satisfying  the  von  Neumann 
stability  criterion.   A  typical  geometry  of  the  problem  is  shown  in 
Figure  2. 
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Boundary  Layer  Equations 

The  Prandtl  boundary  layer  equations  used  in  the  analysis  of 
two-dimensional  mixing  problems  are: 
Global  Continuity 
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where          x  =  u  '^—  (molecular  shear)  (3.6a) 
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J.  =  ~D. (molecular  diffusion)  (3.6b) 
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3y 
3h. 


=  —  T  a. 1-  y  h.J.   (molecular  heat  flux)    (3.6c) 

C  ^      l  3y    H      l  l 
pi         i 

If  the  differentiations  in  Equations  (3.2),  (3.4)  and  (3.5)  are  carried 
out,  the  global  continuity  equation  can  be  extracted  and  the  equations 
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reduce  to  their  more  conventional  form.   The  equations  are  written 
in  the  above  forms  for  the  sake  of  mathematical  expediency  in  the  time- 
averaging  technique  necessary  to  obtain  the  turbulent  forms. 

Inherent  in  the  above  equations  are  the  assumptions  that  the 
body  force  terms  in  the  x- and  y-mcmentum  equations  are  small,  the 
species  are  inert  chemically  (i.e.,  non-reacting  flow)  and  thermal 
diffusion  is  negligible.   From  the  first  of  these  assumptions,  the 
only  information  available  from  the  y-momentum  equation  is  the  result 
that  pressure  is  invariant  in  the  lateral  direction.   Pressure  gradients 
(if  any)  are  restricted  to  the  axial  direction  and  no  further  use 
of  the  y-momentum  equation  can  be  made. 

Further  assumptions  necessary  for  the  present  analysis  are 
summarized  below: 

i.   Flow  outside  the  mixing  zone  is  inviscid,  non- 
conducting and  uniform, 
ii.   Boundary  layers  on  the  confining  walls  will  be  neglected. 
No  correlation  will  be  attempted  for  regions  downstream 
of  the  point  of  interaction  of  the  wall  with  the  mixing 
zone, 
iii.   The  mixing  region  of  interest  is  under  isobaric  conditions; 
that  is,  in  addition  to  oP/3y  =0,  it  is  assumed  that 
9P/3x  =  0.   Although  some  streamwise  pressure  changes  are 
expected  since  the  experimental  configuration  under  considera- 
tion is  one  of  confined  supersonic  flow,  these  pressure 
changes  are  expected  to  be  small  (i.e.,  5-8  percent  of  the 
exit  plane  value)  since  the  region  of  interest  is  the 
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initial  mixing  region  (i.e.,  5-10  inches  downstream  of  the 
nozzle  exit)  where  there  is  no  interaction  of  the  confining 
walls  with  the  mixing  zone.   Data  from  the  experimental 
investigation  justify   this  assumption  as  does  other 
literature  [1] . 

Turbulent  Equations 

In  the  analysis  of  turbulent  flows,  it  is  customary  to  assume 
that  the  instantaneous  value  of  each  property  is  the  sum  of  the  mean 
value  which  varies  with  a  time-mean  average  value  and  a  fluctuating 
component  which  is  a  function  of  time  [3].   In  addition  to  fluctuations 
of  velocity,  density,  temperature  and  mass  fractions,  Van  Driest  [60] 
hypothesized  that  there  are  fluctuations  of  mass  f low  (pu) ,  (pv) 
regarded  as  a  single  property.   Hence,  the  following  relation  is 
defined  for  any  property  $ : 

<j>  =  J  +   <j>'  (3.7) 

where  "bar"  quantities  are  time-mean-averages  defined  by 


$  =  ^-  \  *dt  (3,8) 


and  the  "prime"  quantities  are  fluctuating  components.   It  is  readily 
seen  from  Equation  (3.8)  that  the  time-mean  value  of  any  linear 


fluctuating  quantity  and  its  derivatives  vanislesji.e. ,  u' ,  v1,  3u'/3y, 
etc.  are  all  zero. 

Since  turbulent  flow  is  unsteady  in  nature,  it  is  somewhat  of 
a  contradiction  in  terms  to  speak  of  a  "steady  turbulent  flow",  but 
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the  term  has  usually  been  used  to  denote  a  turbulent  flow  which  is 
steady  in  the  mean,  i.e.,  quasi-steady.   Moreover,  at  each  point  the 
fluid  properties  and  velocity  may  be  observed  to  fluctuate  wildly,  but, 
when  averaged  over  time  periods  comprising  many  cyclic  fluctuations, 
the  time-mean  properties  are  constant  with  respect  to  time. 

Thus,  a  state  of  quasi-steady  turbulent  flow  is  assumed;  the 
instantaneous  values  of  the  quantities  are  replaced  by  the  relation 
in  Equation  (3.7)  and  the  flow  is  considered  in  the  mean  by  time 
averaging  as  defined  in  Equation  (3.8). 

Global  Continuity 

!^  +  !^=0  (3.9) 

dX      dy 
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3x      oy   3y   oy  3x 

This  equation  differs  from  the  laminar  counterpart  by  the  last  two 
terms.   The  last  term  is  the  derivative  of  the  turbulent  normal  stress 
and  is  usually  neglected  by  assuming  that  boundary  layer  approximations 
are  valid  for  cojoint  complex  perturbation  quantities,  i.e., 

kWl-h7^  <3-10) 

The  second  term  on  the  right-hand  side  is  the  Reynolds  (or  apparent) 
stress  which  cannot  be  neglected.   Thus,  a  turbulent  stress  is  defined 
as 


r  "  -  (pv) 'u' 


and  the  momentum  equation  takes  the  form: 
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PU3^+PV    37=   3?[T   +   Tt]  (3'n) 
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pu  - +  pv  t — -  =  -7 t—  (pv)  'a!  -  -r—  (pu)  'a! 
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With  the  approximation  of  Equation  (3.10)  and  defining  turbulent  mass 
diffusion  as 

Jlt  =  -  T7vT^T 

the  specie  equation  takes  the  same  form  as  the  momentum  equation: 

pu   *  +  pv_A.|-  [j.  +  j   ]  (3.12) 

dx       3y    dy  L  i    it 

Total  Energy 


pu  ^+  pv  f£  =  •-  r-^  -  ^  +  f-  [ux  +  u't']  -  -jr  (pv)'H*  -  f-  (pv)'H' 
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Although  the  last  term  can  be  neglected  with  the  aid  of  Equation  (3.10), 
the  present  form  of  the  equation  does  not  permit  a  comparison  between 
the  molecular  and  turbulent  transport  properties.   To  overcome  this 
difficulty,  some  assumptions  have  to  be  made.   From  Equation  (3.6a), 
the  third  term  on  the  right  may  be  written  as 

13   -  3    -2   ""72- 
2^157  IU  +U   ] 

,2 
Along  the  lines  of  Van  Driest's  [60]  work,  it  is  expected  that  u   is 

-2        -2 
small  compared  to  u  ;  since  u  is  assumed  to  be  of  the  same  order  as  h, 

2 

u1   can  be  dropped  from  the  equation. 

Nov  an  equivalent  expression  for  the  perturbation  term,  II1  , 
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will  be  formulated.   From  the  definition  of  total  enthalpy  and  assuming 

2     2 
v  <<  u  ,  which  is  a  valid  boundary  layer  assumption: 

2  2 
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Noting  that   £  a.h.  =  h  -  J  a!h!   and 
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It  is  also  noted  that  Equation  (3.13)  satisfies  the  condition  H'  =  0. 
When  Equation  (3.13)  is  multiplied  by  (pv) ' ,  time  averaged  and  the 
third  order  correlations  neglected,  the  result  is 

-  (cu)'H'  =  -  I   h.(ov)'a!  -  u  (pv)'u'  -  [  a.  (pv)'h!   (3.14) 
i  i 

The  first  and  second  terms  on  the  right  may  be  recognized  as  containing 
the  turbulent  diffusion  and  shear  terms,  respectively.   The  third  term 
is  similar  to  q^  in  Equation  (3.6c)  and  is  defined  as  the  turbulent 
"conduction"  term: 


ict--I  a.  (Pv)'h!_ 


Equation  (3.14)  novj  has  terms  analogous  to  the  molecular  terms  and  the 
total  energy  equation  is  written  as 
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—  3  H  ,  —  3  H   9   /  —     ~\.°   /—     —\,9   —  / —   ,  —  \ 

pu  H  +  pv  37  =  ^7  (qct  -  q)  +  a7  (qdt  -  V  +  ^7  u(Tt  +  T) 

(3.15) 
It  is  now  assumed  that  the  boundary  layer  is  fully  turbulent 
and,  thus,  molecular  transports  are  negligible,  i.e., 
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The  turbulent  shear  is  related  to  the  mean  flow  variables 
following  Boussinesq  [61] 

3u 


T  -  -  (pv)'u'  =  e  ~  (3.16) 

t  m  3y 

where  e   is  the  turbulent  momentum  transfer  coefficient  and  has  the 
m 

same  units  as  the  molecular  dynamic  viscosity.   Similarly,  the 
turbulent  diffusion  and  conduction  terms  are  related  to  the  mean  flow 
variables : 

3a. 
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The  relation  indicated  by  Equation  (3.17)  implies  that  all  mass 

transfer  coefficients,  e,.,  of  the  various  species  are  equal.   Woolley 
di 

[17]  argued  that  if  two  streams,  each  composed  of  multiple  chemical 
species  in  homogeneous  mixtures  were  considered,  then,  upon  exposing 


28 


the  two  streams  to  each  other,  all  gradients  of  concentration  differences 
for  species  between  the  streams  would  be  identical.   If  their  mass 
transport  coefficients  were  also  equal,  they  would  diffuse  through  the 
mixing  zone  at  the  same  rate.   The  relative  concentrations  of  the  species 
from  a  given  stream  would,  then,  remain  unchanged  at  any  position  in 
the  mixing  zone.   However,  this  was  equivalent  to  each  stream  behaving 
as  a  single  species.  Thus,  under  the  present  assumption,  each  homogeneous 
stream,  no  matter  what  its  detailed  composition,  may  be  treated  as  a 
single  species  in  the  mixing  study.   Therefore,  it  is  only  necessary  to 
treat  the  mixing  of  two  dissimilar  gases,  each  having  the  average 
chemical  and  thermodynamic  properties  of  their  respective  mixtures. 
This  is  particularly  true  for  the  primary  stream  of  air  which  is 
usually  considered  as  a  single  species [62]. 

Also  of  interest  are  some  dimensionless  turbulent  quantities 
which  measure  the  relative  rates  of  different  transport  mechanisms. 
The  turbulent  Prandtl  number  is  a  measure  of  the  relative  rates  of 
transport  of  momentum  and  energy;  the  turbulent  Lewis  number  is  a 
measure  of  the  relative  rates  of  transport  of  mass  and  energy. 

These  quantities  together  with  the  turbulent  Schmidt  number 
are  defined  as 
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Since  all  equations  are  time  averaged  and  all  transports  are 
turbulent,  the  bar  notation  and  the  subscript  "t"  will  be  dropped 
from  here  on.   With  the  definitions  in  Equations  (3.16)  through  (3.19) 
and  after  some  rearrangement,  the  governing  equations  take  the  final 
form: 
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The  von  Mises  Transformation 

The  solution  of  Equations  (3.20)  through  (3.23)  provides  the 
details  of  the  flow  field  including  the  velocity,  species  and  enthalpy 
(thus  temperature)  fields.   The  global  continuity  equation,  Equation 
(3.20),  can  be  eliminated  from  the  system  of  differential  equations  by 
introducing  the  von  Mises  coordinates  as  the  independent  variables. 
The  transformation  (>:,y)  ■>  (x,^(x,y))  is  defined  according  to  the 
relations: 
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The  derivatives  in  the  physical  coordinates  are  mapped  onto  the 
von  Mises  plane  via 
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Substitution  of  the  above  relations  into  the  system  of  equations 
completes  the  formulation  of  the  problem  in  the  von  Mises  plane. 
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■r—   DUG       -TT~ 

Sc          m   dv 

Total   Energy 

9H  m   3__ 

dX          cv 

m 

1      3H        Pr   -   1   3u2/2 
Pr   difi            Pr          diji 

+  ii 

3r        y      : 

9a. 

i 

L  9^ 

(3.26) 


(3.27) 


(3.28) 


The  physical  y-coordinate  is  obtained  by  the  inverse  transformation: 
* 


dl 


(3.29) 


and    the   transverse   component   of   velocity,   v,    is   given  by; 


p     dX 


(3.30) 
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Boundary  Conditions 

The  governing  equations  exhibit  parabolic  characteristics 
and  thus  require  initial  conditions  at  some  x  -  x  and  boundary- 
conditions  at  ty  ~   0  and  if;  =  °°. 

In  any  real  experimental  situation,  there  is  an  inevitable 
accumulation  of  boundary  layer  on  the  jet  dividing  boundary.   In 
order  to  avoid  making  an  error  in  initial  conditions  by  assuming 
either  a  step  profile  or  a  computed  boundary  layer,  the  calculation 
of  the  mixing  region  is  started  at  a  position  downstream  of  the  mixing 
interface  where  measured  data  are  available.   Hence,  the  initial 
conditions  may  be  expressed  as: 

@x  =  x    ;  Q    <_  ty   <_  <*> 

*       * 

u(x  ,l|0  =  u  d>) 

(3.31) 
H(x*,ijO  =  H*(U>) 

a.(x*,^)  =  a*W 

Since  the  wall  boundary  layer  is  assumed  to  be  negligible, 
constant  flow  conditions  equal  to  the  secondary  stream  conditions  are 
assumed  at.  the  wall,  i.e.,  slip  condition.   It  is  observed  that  zero 
axial  velocities  are  not  permissible  at  any  location  in  the  flow  field 
due  to  the  inverse  transformation  of  Equation  (3.29).   Then,  the 
boundary  conditions  become: 

8  y  =  0  -  tj,  =  0  ;  x*  <  x  <  xmax 
u(x,0)  =  u  =  constant 
H(x,0)  ■  H  -  constant 
u0  (x,0)  =  0  ;  Oj,  (x,0)  -  0  ;  aA)He(x,0)  =  1  ' 


(3.32) 
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(3.33) 


@  Y 


u(x,^)  =  u  =  constant 
P 

H(x,°0  =  H  =  constant 
P 

yx,-)  =  0.23  ;  a^Cx.co)  -  0.77  ;  ^(x,-) 


The  condition  of  y  =  «  (ifi  =  «)  are  all  regions  beyond  the  upper 
boundary  of  the  mixing  zone. 

Finite  Difference  Equation': 

A  forward  marching  all-explicit  numerical  method  was  used  in 
this  analysis.   Accordingly,  for  the  transverse  derivatives,  central 
differences  are  used  in  the  interior  grid  points.   Forward  differences 
are  used  for  the  longitudinal  derivatives  everywhere.   Figure  3  shows 
a  generic  point  (n+l,tn)  in  the  (x,">»  grid  network,  for  which  the  solution 
is  obtained  by  using  the  following  explicit  difference  relations  [63] 


(3.34) 


(3.35) 

^rt,m-H/2   n,m-fl    n,m ~n,m-l/2   n,m    n,m-l 

(3.36) 


3F   Fn+l,m 

F 
n,m 

Sx        Ax 

3F   Fn,mfl 

F     . 
n,m-l 

3y       2;., 

':  V 


5^ 

4  H 


where  ?  =   puc  (3.37) 

m 

^n,m±l/2  "  1   (S>,m  +  ^.m+P  (3'38) 

Then,  the  non-linear  parabolic  equation  may  be  approximated 
by  the  finite  difference  equation 
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Fn+l,m       Fn,m  +  y[ ?n,m+l/2Fn,m+l  '    (?n,m+l/2  +   ^n,m-l/2)Fn,m 


+  S  inF  .]    +  (3.39) 

n,m-l/2  n,m-l 

where  y  -  Ax/(A',j/)   and  £  is  the  error  introduced  by  the  finite  difference 
approximation  of  the  differential  equation  and  may  he  defined  as  [63]: 


£  =  k1[e(Ax)]  +  k2[e(£^)2]  (3.40) 

Equation  (3.40)  implies  that  Ax  and  Aip  have  to  be  sufficiently  small 
for  the  difference  scheme  to  be  accurate.   This  point  will  be 
discussed  further  with  respect  to  consistency  and  stability  of  the 
scheme. 

It  should  be  recognized  that,  although  the  partial  differential 
equations  are  non-linear,  the  present  explicit  difference  formulation 
results  in  a  locally  linear  system.   Inherent  in  this  result  is  the 
assumption  that  the  solution  is  "fairly  smooth"  and  the  quantity,  t, , 
is  a  "slowly  varying"  function.   These  two  assumptions  imply  that 
discontinuities  such  as  shock  waves  cannot  be  present  in  the  flow 
field,  and  thus  pose  no  major  restriction  on  the  solution  since  no 
such  phenomena  are  considered. 

The  conservation  equations  for  the  interior  grid  points 
(i.e.,  m  -■£   0)  in  difference  form  are: 
Momentum 


Vfl.m  ""  Un,m  +  Y  [  Cn,m+l/2Un,m+l    (Cn,mfl/2  +  Cn,m-l/2)un  ,m 


+  K         1/9u    J  (3.41) 

n,m-l/2  n,m-l 
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Species 

i''n+l,m        vv*i/n,ra    '    Sc    L  ~'n,m+l/2  ^"'i'njm+l        x-°n,m+l/2 


£  ,  /0)(cO  +   £  .  /0(a.)  J  (3.A2) 

n,m-l/2        1  n,m  n,m-l/2      1   n,m-l 


Hn+l,m        Hn,m  +  Pr    [Sn,m+l/2Hn,nrf  1        (Sn,m+l/2  +   5n,ir,-l/2)Hn,i 

+   E  HI     I    ^Pr    -   1>    rr  „2 

Si,m-l/2   n,m-lJ  2Pr  L^n,ro+l/2  n,m+l 


2  2 

(Cn,m+l/2  +   Cn,m-l/2)Un,m  +   ?n,m-l/2Un,m-l] 


^TT21  t(»i)n,n,+l/2'»i>n,^l  "    K»i>„,«/! 


l  n,m-l/Z        l  n,m  in, m- 1/2      i  n,m-l 

(3. A3) 

Boundary  Conditions  for  the  Difference  Equation 

The  initial  and  boundary  conditions  for  the  difference 
equations  are  similar  to  their  counterparts  for  the  differential 
equation  and  are  input  in  equal  intervals  of  AtJ/. 

@x  =  x   ;  0  <_  ij;  <_  NA^j 

u    =  u  (AiJ/) 
°'m    m  (3.44) 

i  o,m     i     m 
An  interesting  characteristic  of  the  general  difference  equation, 
Equation  (3.39),  is  observed  when  the  boundary  condition  at  iJj  =  0,  i.e., 
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m  =  0,  is  applied;  that  is,  the  terms  F   ,  and  ?n  ,  .^   are  undefined. 

The  conventional  method  employed  to  circumvent  this  difficulty  is  to 

define  an  "artificial"  boundary  condition  such  that  the  value  of 

F   ,  is  equal  to  F   . .   This  implies  that  the  gradient  of  the  quantity 
n,-l  n ,  1 

F  is  zero  at  the  boundary  and,  in  the  physical  sense,  is  equivalent  to 
a  condition  of  symmetry  at  the  axis  of  a  jet.   Since  an  assumption  of 
constant  velocity  at  the  wall  has  already  been  made  for  the  physical 
boundary  condition,  the  application  of  the  above  principle  will  not 
introduce  any  new  assumptions  into  the  analysis.   With  the  application 
of  the  condition  of  symmetry,  the  difference  equations  at  the  wall 
become: 

Momentum: 

u  ,,    -  u    +  2yC    (u  .  -  u   )  (3.45) 

n+1 ,o    n,o     '  n ,  o  n,l    n,o 

Species 

(a.)  ^   -(a.)    +|^C   [(a.)   ,-Ca.)   J  (3.46) 

i  n+ljO     l  n,o   Sc  n,o   l  n,l     i  n,o 

Energy 

Hn+l.o  =  Hn,o  +  £  ^n,o<Hn,l  "  \^   +  ^^   ^n,o(un,l 


2y(Le  -  1)      £  (   >    [(   }    .  (   s 

'  Pr       n,o  L        i  n,o   i  n,l     l  i 


(3.47) 


Together  with  the  boundary  condition: 


H  =   H      =   constant  (3.48) 

n,o  s 

<°02>n.o  =  °    ;    (aN2}n,o  =  °   ;    <aA,He>n,o  =  X 
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The  boundary  condition  at  t|>  =  »  is  one  of  "floating"  type;  that  is, 

calculations  proceed  until  the  condition  F    .,  ■  F    is  satisfied. 

n,m+l    n,m 

Analysis  of  the  Difference  Scheme 

Once  a  finite  difference  scheme  is  set  up,  it  must  satisfy 
three  conditions  [64] : 

i.   The  difference  equation  must  be  "consistent"  with 

the  differential  equation.   That  is,  the  error  involved 
in  approximating  a  differential  equation  by  a  difference 
equation  must  vanish  in  the  limit  Ax,  A\J>  ■*  0. 
ii.   The  difference  scheme  must  be  "convergent".   That  is,  the 
solution  to  the  difference  equation  must  approach  the 
solution  of  the  differential  equation  in  the  limit. 
iii.   The  marching  rate  must  be  "stable".   That  is,  the  grid 
size  should  be  such  that  the  solution  does  not  become 
unbounded  anywhere.   A  relation  between  Ax  and  A^  must 
be  found  to  insure  this  condition. 
The  consistency  of  the  present  method  may  be  shown  by 
observing  that  the  quantity  <S,  defined  in  Equation  (3.40),  approaches 
zero  unconditionally  as  (Ax,  A4-)  approach  zero.   The  proof  of  conver- 
gence for  quasi-linear  difference  equations  is  lengthy  and  is  not 
presented  here;  however,  the  numerical  scheme  employed  can  be  shown 
to  ba  convergent  with  the  method  of  Strang  [65]. 

One  of  the  important  characteristics  of  explicit  finite 
difference  schemes  is  the  stability  of  the  solution.   For  a  system  of 
linear  constant  coefficient  equations,  it  can  be  shown  that  the 
von  Neumann  stability  criterion  is  the  necessary  and  sufficient  condition 
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for  stability  and  convergence  [63].   In  variable-coefficient  cases  it 
is  necessary  for  stability  that  the  von  Neumann  conditions  (derived  as 
though  the  coefficients  were  constant)  be  satisfied  at  every  point 
in  the  grid  network.   This  is  based  on  the  observation  that  when 
instability  occurs  in  practice,  it  often  appears  as  a  local  disturbance 
in  a  region  where  the  von  Neumann  condition  is  violated.   John  [in  63]  h; 
shown  that  for  a  general  class  of  explicit  difference  equations  for 
quasi-linear  parabolic  partial  differential  equations,  the  von  Neumann 
condition  is  necessary  for  stability  and  that  a  slightly  modified  form 
of  the  von  Neumann  condition  is  sufficient. 

On  this  basis  the  von  Neumann  stability  critericnwas  derived, 
and  this  critericnwas  checked  at  each  point  in  the  grid  network  since 
the  system  of  governing  equations  are  quasi-linear.   An  outline  of  the 
von  Neumann  analysis  for  the  difference  equations  used  is  given  below. 

A  linear,  constant  coefficient  parabolic  partial  differential 
equation  with  periodic  boundary  conditions  has  a  solution  of  the  form 

u(x,y)  =  I   \(x)eiky 
k  k 

Assuming  that  in  the  limit  the  difference  equation  has  the  same  type 
of  solution,  let 

*    k 
Equation  (3.49)  is  substituted  into  the  general  difference  equation, 
Equation  (3.39),  and  since  it  is  also  assumed  that  all  Fourier 
coefficients  decay  exponentially,  the  k   term  is  examined: 

\=1~   Y(?n,m+l/2  +  f'n,m-l/2)(1  "  C0S  kA*>  +  ^n,m+l/2 

-  5    ,  /9)sin  kAy  (3.50) 

n  ,m-i/  /. 
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For  F    to  be  bounded,  the  von  Neumann  criterion  is 
n,m 

M  <  i 

After  some  manipulation,  the  restriction  on  the  grid  size  relation. 
Y,  is  obtained 

gn,m+l/2  +  Cn,m-l/2 


2  2 

5         ■  i  in  +  €  t  /o   -   2C         ,,  /0C  ,  /0    cos   kAy 

n,ra+l/2    n,m-l/2     n,m+l/2  n,m-l/2 


and  the  most  restrictive  condition  occurs  when  the  quantity  on  the 
right  is  a  minimum,  which  is  the  case  when  cos (kAy)  =  -1.   Then,  the 
final  form  of  the  stability  criterion  which  governs  the  interior 
grid  size  is 

Ax  1 


(A*)     Si,m+l/2   Cn,m-l/2 


U   =   pue  )      (3.51) 


It  is  further  observed  that  if  ;  is  a  constant,  C,  throughout 

2 

the  domain  of  interest,  Equation  (3.51)  reduces  to  the  form  CAt/(Ax) 

1/2  which  is  the  well-established  stability  condition  for  the  linear 
"heat  equation"  with  the  difference  scheme  of  Equation  (3.39). 

Since  Equation  (3.51)  is  valid  only  at  the  interior  grid 
points,  a  stability  condition  for  the  boundary  is  derived  in  the 
same  manner  with  the  result 


(3.52) 


(AijO      n,o 

The  above  stability  conditions  were  derived  for  the  general 
difference  equation.   V.Tien  applied  to  the  governing  system  of  equations, 
they  are  seen  to  be  different  by  a  constant  K  where 


:j/j 


K  =  1  for  the  momentum  equation 

K  =  1/Sc       for  the  species equati on 
K  =  1/Pr       for  the  energy  equation 

It  was  found, during  the  analysis  of  truncation  error  of 
linear  equations,  that  when  the  stability  criterion  is  multiplied  by 
a  factor  of  1/3,  a  higher  order  of  accuracy  could  be  obtained,  i.e., 

c  2  4 

C  =  K  [6 (Ax)  ]  +  K  [e(A'^)  ].   This  is  not  strictly  correct  for  quasi- 
linear  equations,  but  if  £  is  a  "slowly  varying"  function,  considerable 
improvement  in  accuracy  can  still  be  obtained  by  making  use  of  the 
above  result.   Therefore,  with  the  above  modifications  the  final  form 
of  the  stability  conditions  are 

iiK^U  Jft2e    1/2> 

n,o  n,m+l/2    n,m-l/2 

This  results  in  six  conditions,  the  most  restrictive  of  which  is 
utilized,  depending  on  the  magnitude  of  the  dimensionless  quantities 
Pr,  Sc  and  Le. 

Thus,  Equations  (3.41)  through  (3.48)  together  with  the  stability 
condition,  Equation  (3.53),  constitute  the  numerical  solution  to  the 
flow  field.   The  only  remaining  point  is  the  formulation  of  the 
quantity  e  ,  and  this  will  be  discussed  in  Chapter  IV. 


CHAPTER  IV 
TURBULENT  MOMENTUM  AND  MASS  TRANSPORT 

General  Considerations 

Although  the  first  portion  of  this  chapter  could  have  been 
included  in  the  literature  survey,  it  is  felt  that  the  presentation 
of  the  analysis  would  be  more  continuous  if  it  were  included  in  this 
section. 

To  this  point  the  consideration  of  the  analytical  treatment 
of  the  turbulent  mixing  problem  has  dealt  primarily  with  the  question 
of  solving  the  questions  of  motion.   While  being  necessary,  this  is 
not  the  area  of  greatest  difficulty.   The  major  difficulty  associated 
with  these  problems  is  the  mathematical  representation  of  turbulent 
transport  processes.   The  present  understanding  of  the  turbulence 
phenomenon  is  such  that  turbulent  processes  within  a  shear  flow  cannot 
be  treated  locally.   Rather,  the  most  that  can  be  expected  is  some 
prediction  of  the  "mean"  flow  properties. 

There  are  essentially  two  distinct  approaches  previously  taken 
in  the  analysis  of  mixing  problems.   The  more  recent  approach  [66,67,68] 
is  to  introduce  additional  conservation  equations  which  describe  the 
Reynolds  stress.   This  approach  is  appealing  since  it  guarantees 
conservation  of  the  turbulent  quantities.   Unfortunately,  the  resulting 
equations  contain  second  and  higher  order  correlations  and,  thus,  to 
apply  this  approach,  empirical  relations  are  required  for  the  third 
and  higher  order  correlations.   Furthermore,  due  to  its  complexity, 
this  approach  has  not  been  applied  to  compressible  free  shear  layer 


AO 
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flows.   Therefore,  solutions  of  practical  problems  of  current  interest 
have  been  attained  only  by  employing  the  more  commonly  used  approach 
of  utilizing  an  eddy  viscosity  model  along  with  assumed  constant  Sc, 
Pr  and  hence  Le.   The  eddy  viscosity  models  that  have  been  proposed 
for  free  mixing  flows  are  summarized  in  Table  I. 

Classical  Eddy  Viscosity  Models 

The  famous  mixing  length  theory  for  turbulent  shear  was 
formulated  by  Prandtl  [7]  who  hypothesized  that  the  mean  value  of  the 
fluctuating  velocity  component  in  a  turbulent  flow  field  is  equal 
to  the  product  of  the  local  mean  velocity  gradient  and  a  characteristic 
mixing  length,  1.      The  quantity  I   is  defined  as  a  distance  in  the  flow 
field  such  that  a  fluid  element  conserves  its  longitudinal  velocity 
as  it  moves  across  this  distance.   In  the  case  of  free  mixing,  the 
mixing  length  is  assumed  to  be  constant  across  the  mixing  layer  and 
also  assumed  to  be  proportional  to  the  local  width  of  the  mixing  zone. 
Thus,  Prandtl 's  mixing  length  theory  gives  the  following  relation 
for  the  Reynolds  stress  in  the  longitudinal  direction 

-; r-| — r   „   2,2   3u    3u  ,.  1N 

t  =  -  (pv)'u'   Ub   —   r-  (4.1) 

t  |  dy  I  3y 

where  c  is  an  experimentally  determined  constant  and  b  is  the  width 
of  the  mixing  layer.   Using  the  concept  of  eddy  viscosity  in  Equation 

(3.16) 

2,2   3u  I  //on 

e   =  c  b   —  (4.2) 

m       |  3y  I 

Based  on  a  similar  mixing  length  concept,  Taylor  [in  1]  derived  a 
vorticity  transport  theory  where  the  vorticity  of  the  fluid  element 
is  assumed  to  be  conserved  across  the  mixing  length.   In  both  of  these 
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mixing  length  concepts,  the  eddy  length  scale  was  assumed  to  be  much 
smaller  than  the  local  width  of  the  mixing  layer.   The  complexity  of 
Taylor's  model  for  axi--symmetric  flows  has  prevented  its  utilization; 
and  for  the  case  of  two-dimensional  flows,  except  for  a  numerical 
constant,  the  vorticity  theory  results  in  the  same  expression  for 
turbulent  shear  as  obtained  from  Prandtl's  mixing  length  theory. 

Prandtl  [10]  later  proposed  another  model  for  the  turbulent 
eddy  viscosity  based  on  the  hypothesis  that  the  eddy  scale  was  of  the 
same  order  as  the  width  of  the  mixing  layer.   This  model  was  based  on 
the  assumption  that  the  eddy  viscosity  is  related  to  the  local  mean 
velocity  gradient  and  was  expressed  as 

e  =  cb(U    -  U  .  )  (4.3) 

m       max    mm 

where  c  is  an  empirical  constant  and  U    and  U  .   are  the  maximum  and 

max      mm 

minimum  longitudinal  velocities,  respectively.   Prandtl's  second  model 
predicts  a  constant  eddy  viscosity  across  the  mixing  layer  since  the 
width,  b,  is  not  a  function  of  the  lateral  coordinate. 

Equation  (4.3)  has  been  widely  applied  to  a  variety  of  free 
mixing  problems  due  to  its  mathematical  simplicity  and  the  results  it 
yields  agree  satisfactorily  with  experimental  data  for  several  flow 
configurations.   It,  therefore,  forms  the  basis  for  most  eddy 
viscosity  models  existing  in  the  literature,  and   is  used  for  a 
particular  flow  field  by  appropriately  including  the  effects  that  may 
be  of  significant  interest  in  that  case.   However,  it  is  noted  that 
the  model  fails  completely  when  the  velocities  of  the  two  streams  are 
equal.   Hence,  it  predicts  that  two  streams  of  equal  velocities  flow 
along  as  segregated  without  turbulent  mixing.   This  implication  has 
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been  shown  to  b   incorrect  by  the  experimental  results  of  references 
[5]  and  [35]. 

The  various  eddy  viscosity  models  available  today  are  the 
result  of  attempts  by  several  investigators  to  include  the  cases  of 
equal  velocities  within  the  framework  of  Prandtl's  original  hypothesis 
for  free  turbulent  flow.   Ferri  et  al.  [34]  suggested  the  following 
model  by  simply  extending  the  second  Prandtl  model  to  describe  flows 
with  density  gradients 

m      1//     max       mm 

where  r  ,„  is  the  half-radius.   This  model  has  yielded  predictions  of 

unreliable  accuracy  for  axi-symmetric  flows  [35,58]  but  when  applied- 

to  the  planar  case,  good  predictions  were  achieved  [36].   However,  the 

Ferri  model  fails  when  the  mass  flux  of  each  stream  is  equal. 

In  order  to  circumvent  such  irregularities,  Alpinieri  [35] 

considered  the  eddy  viscosity  to  be  proportional  to  the  sum  of  the 

mass  flux  and  the  momentum  flux  and  for  an  axi-symmetric  jet,  proposed 

the  relation 

2- 

(4.5) 


(pSn>C.L. 


(pu)j 


2 


where  (pe  )  is  the  "dynamic"  turbulent  viscosity  at  the  center-line 
m 

of  the  coaxial  jet  with  the  subscripts  e  and  j  referring  to  the 
properties  of  the  external  stream  and  the  inner  jet,  respectively. 
The  Alpinieri  model  is  contrary  in  form  to  any  other  model  and  is 
viewed  as  essentially  empirical,  qualitatively  as  well  as  quantitatively 
[69]. 
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Density  differences  may  arise  within  the  mixing  region  either 
due  to  compressibility  effects,  as  in  the  case  of  heated  jets  and 
supersonic  flows,  or  due  to  streams  of  different  composition.   Ting 
and  Libby  [70],  employing  a  Mager  transformation,  postulated  the 
following  relation  between  the  eddy  viscosity  for  constant  density 
mixing  and  that  for  variable  density  axi-symmetric  flows 


I^M'fc- 


where  e   is  the  eddy  viscosity  for  incompressible  flows  and  p   is  a 
reference  density.   As  can  be  seen,  the  above  relation  is  essentially 
a  conversion  of  the  incompressible  eddy  viscosity,  £  ,  to  one  applicable 
for  flows  with  density  variations  either  due  to  compressibility  or 
stratification.   It  should  be  recognized,  however,  that  while  this 
transformation  admits  possible  practical  applications,  no  definite 
form  of  c  or  p   is  suggested  and  the  results  vary  depending  on  the 
forms  of  e  and  p  used.   A  planar  form  of  this  model  was  utilized 
by  Schetz  [36]  with  the  expression 


P2^  : 


I-   > 


r 
*  2     1    p   . 
:  pcl.  hjTj  pT-'1 


(A. 7) 


Donaldson  and  Gray  [20]  attempted  to  account  for  compressibility 
via  modifying  the  Prandtl  model,  Equation  (4.3),  by  employing  an 
empirically  determined  constant  which  varied  with  Mach  number  at  the 
half-radius.   This  resulted  in  the  expression 

e  =  £  (0.66  +  0.34  exp(-3.42  M2))  (4.8) 

m    p 
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where  e  denotes  the  Prandtl  eddy  viscosity  model  of  Equation 
(4.3). 

Schetz  [69]  proposed  a  model  which  was  developed  from  an  exten- 
sion of  Clauser's  model  for  the  wake  region  of  a  turbulent  boundary 
layer  to  free  shear  layers.   The  specific  functional  expression  for 
a  given  flow  problem  was  derived  from  the  general  statement:   "the 
turbulent  viscosity  is  proportional  to  the  mass  flow  defect  (or 
excess)  in  the  mixing  region."  [69,  page  1]   The  model  was  expressed  as 

r  I         | 

pc   =  cp  u   Ml  -  2 dy  (4.9) 

m  e  e  J  !     Peue  I 

Similar  to  the  Ferri  model,  the  above  expression  fails  when  the  mass 
defect  (or  excess)  is  zero.   Also,  the  model  has  been  shown  [31]  to 
fail  for  the  case  of  the  quiescent  jet  since  the  mass  entrainment, 
and  therefore,  mass  defect,  increases  in  the  downstream  direction 
thus  predicting  continually  increasing  eddy  viscosity.   It  has  been 
demonstrated  by  Eggers  [71]  that  for  an  accurate  prediction  of  the 
flow  field  in  quiescent  supersonic  jets,  the  eddy  viscosity  must 
remain  very  nearly  constant. 

There  have  been  other  models  formulated  for  the  turbulent 
transport  mechanism,  but  these  are  either  too  complicated  to  be 
of  practical  use,  i.e.,  von  Karman  model  (Table  I),  or  are  not 
applicable  to  the  present  study,  i.e.,  Zakkay  model.   Each  of  the  eddy 
viscosity  models  presented  in  Equations  (4.2)  through  (4.9)  are 
deficient  in  certain  respects.   Of  course,  some  of  the  deficiencies  are 
relatively  unimportant  and  any  one  of  these  models  is  acceptable 
provided  that  they  are  used  within  the  flow  region  for  which  they  have 
been  verified. 
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Choice  of  Models  for  the  Present  Studv 


It  should  be  mentioned  that  the  eddy  viscosity  models  cited 
have  been  used  for  extensive  correlations  but  only  in  the  similarity 
regions  of  symmetric  jets,  subject  to  specific  flo\v'  conditions  and 
configurations.   It  is  important  to  note  the  distinction  between  the 
similarity  region  of  a  symmetric  jet  and  the  similarity  region  of  the 
half-jet  which  is  the  case  under  study.   The  similarity  region  of  a 
symmetric  jet  is  that  part  of  the  flow  field  far  downstream  of  the 
so-called  potential  core.   The  half-jet  may  be  thought  of  as  a  symmetric 
jet  with  an  infinite  jet  radius  or  height,  resulting  in  a  potential 
core  of  infinite  length.   The  important  point  is  the  fact  that  the 
models  in  Equations  (4.4)  through  (4.9)  have  been  verified  only  in 
regions  past  the  potential  core  of  a  symmetric  jet.   Thus,  no  state- 
ment may  be  made  regarding  their  application  to  other  jet  mixing 
configurations.   In  fact,  this  is  true  for  all  available  expressions 
for  e  because  of  the  lack  of  complete  and  accurate  data  used  in 
studying  any  eddy  viscosity  model.   Hence,  the  results  obtained  in 
this  study  apply  to  both  the  initial  and  fully  developed  regions  of  a 
half-jet  and  the  potential  core  of  a  symmetric  jet. 

The  fact  that  several  models  may  be  used  to  correlate  the  same 
experimental  data  is  shown  by  the  results  obtained  by  Ragsdale  and 
Edwards  [72,73]  in  their  analytical  and  experimental  study  with  air- 
bromine  system.   In  their  analytical  study,  various  expressions  for 
eddy  viscosity  were  compared  on  a  consistent  basis.   It  was  concluded 
that  modifications  of  Prandtl's  second  hypothesis,  that  introduce  mass 
flux  or  momentum  flux  or  both,  produce  expressions  whose  differences 
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are  more  apparent  than  real.   It  was  shown  that  these  various  expressions 
predict  essentially  the  same  eddy  viscosity  as  long  as  they  are  applied 
only  within  the  range  of  conditions  for  which  they  have  been  experimentally 
verified.   It  was  concluded  that  this  was  perhaps  because  the  initial 
turbulence  present  in  the  streams  contributes  significantly  to  the 
mixing  process  and  may  dominate  the  situation  for  nearly  equal  stream 
velocities. 

The  effect  of  free  stream  turbulence  was  also  considered  by 
Hokenson  and  Schetz  [74]  in  their  study  of  turbulent  mixing  with 
pressure  gradients.   The  results  of  their  investigation  demonstrated 
that  the  empirical  constant  in  the  modified  Clauser  model  of 
Equation  (4.9)  must  reflect  the  turbulence  intensity  and  therefore  is 
not  a  universal  constant.   It  was  also  observed  that  this  empirical 
constant  was  approximately  independent  of  the  longitudinal  distance 
from  the  initial  station,  and  that  if  adequate  information  (i.e., 
turbulence  intensity)  is  known  at  the  initial  station,  a  numerical 
evaluation  of  the  flow  field  could  reasonably  be  assumed  within  the 
framework  of  the  generalized  Clauser  eddy  viscosity. 

To  correlate  the  data  obtained  in  the  experimental  phase  of  the 
present  investigation,  four  models  of  the  turbulent  momentum  transport 
mechanism  were  chosen;  namely,  Prandtls  mixing  length  hypothesis 
(Equation  (4.2)),  Ferri's  differential  mass  flux  model  (Equation  (4.4)), 
Alpinieri's  momentum  flux  model  (Equation  (4.5))  and  Schetz' s 
extension  of  the  Clauser  model  (Equation  (4.9)).   As  stated  earlier 
these  models  have  been  shown,  within  certain  limitations,  to 
correlate  data  in  the  "far  field"  or  similarity  region  of  symmetric 
jets.   Since  the  region  of  interest  for  the  present  analysis  is  the 
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initial  mixing  region,  the  primary  objective  was  to  examine  the 
validity  of  these  models  in  the  so-called  "near  field".   It  is  also 
of  passing  interest  to  note  that,  to  the  author's  knowledge,  the 
Alpinieri  model  has  never  been  applied  to  the  case  of  two-dimensional 
mixing  problems  in  any  region  of  flow.   Minor  modifications  necessary 
to  compare  the  models  on  a  consistent  basis  within  the  framework 
of  the  present  analysis,  is  discussed  in  the  next  section. 

Modification  and  Correlation  of  Models 


The  form  of  the  governing  equations,  Equations  (3.26)  through 

(3.30),  require  that  the  turbulent  transport  mechanism  be  specified  in 

the  form  (or  units)  of  "dynamic"  turbulent  viscosity;  i.e.,  as  the 

counterpart  of  u  in  the  laminar  case.   To  avoid  confusion  in  terms  of 

notation,  (u  )„,  (u  )_,  (u  ),  and  (u  )_  will  be  used  to  denote  the 
tr     tr     tA        tb 

Prandtl,  Ferri,  Alpinieri  and  Schetz  models,  respectively. 

To  be  consistent  with  the  governing  equations,  the  Prandtl 
model  is  expressed  in  terms  of  the  von  Mises  variables 

(VP  =  ^  =  cib2p2u|  1^1  (A-10) 

and  with  the  aid  of  Equation  (3.35)  is  put  in  finite  difference  form 

GO-  =  c?b2  lD2   u  (4.11) 

t  P    1  n-1  n,m  n,m      2Aip 

In  the  Ferri  and  Alpinieri  models,  the  eddy  viscosity  is 
expressed  in  terms  of  the  half-radius  based  on  either  the  velocity  or 
the  mass  flux.   The  half -radius  is  defined  as  the  distance  from  the 
axis  of  symmetry  at  which  the  axial  velocity  (or  mass  flux)  is  equal  to 
the  average  of  the  maximum  and  minimum  velocities  (or  mass  flux),  i.e., 
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r  at  which  U  =  0.5 (U    +  U  .  ) 

max    mm 

r  at  which  pU  =  0.5((pU)    +  (pU)  .  ) 

max       nun 


Although  this  has  some  physical  meaning  in  the  similarity  region  of  a 
symmetric  jet,  it  has  no  meaning  in  the  initial  region  of  a  half-jet. 
Therefore,  rather  than  a  width  such  as  the  half-radius  being  used,  the 
actual  width  of  the  mixing  region  will  be  utilized.   Since  there  is 
always  an  uncertainty  as  to  the  physical  boundaries  of  the  mixing 
region  in  any  experimental  study,  the  width  will  be  defined  as  the 
distance  in  which  the  velocities  are  within  5  percent  of  the  free 
stream  values,  i.e.,  referring  to  Figure  2,  if 

b  -  Yl  "  Y2 

then      Y.  =  Y   at  which   I U  -  U  I /U  =0.05 

1  '     P1   P 

Y  =  Y   at  which   lu  -  U  l/U  =0.05 

2  '     s '   s 


(4.12) 


In  the  Alpinieri  model  the  centerline  velocity  is  replaced  by  the 
secondary  stream, which  is  analogous  to  the  condition  in  the  potential 
core  of  the  symmetric  jet. 

With  the  above  modifications,  the  Ferri  and  Alpinieri  models 
take  the  form 


(ujv   =  c„b[(pu)    -  (pu)  .  ] 
t  F    I  max       mm 

2-1 

(Oa  =  c3b(pu)s 


P    P  u 

-JL  +  _E_JL 

P       2 

s   p  u 

Hp  s- 


(4.13) 


(4.14) 


where  the  subscripts  p  and  s  refer  to  primary  and  secondary  stream 
conditions,  respectively. 

The  Schetz  model  when  expressed  in  terms  of  the  von  Mises 


variables  reduces  to 
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P  u 

1£_E 


(Vs  =  CA      "^f-1  d*  (4'15) 


It  was  necessary  to  have  a  common  frame  of  reference  to  be 
able  to  compare  the  different  models  on  a  consistent  basis.   The 
mixing  zone  width,  defined  in  Equation  (A. 12)  was  chosen  to  be  this 
common  reference.   The  coefficients  of  each  model  were  varied  until 
the  predicted  growth  of  the  mixing  zone  matched  the  experimental  data 
for  each  of  the  four  configurations  of  Table  II.   With  each  model 
predicting  the  same  growth  rate,  the  theoretical  and  experimental 
velocity  profiles  were  compared  as  to  the  "goodness  of  fit"  of  each 
model. 

It  was  also  expected  that  the  coefficients  of  a  given  model 
would  vary  from  configuration  to  configuration  and  thus  be  a  function 
of  the  flow  field  as  was  concluded  by  Hokenson  and  Schetz  [74], 
However,  rather  than  use  the  turbulence  intensity  at  the  initial 
station,  an  attempt  was  made  to  correlate  the  coefficients  (C   through 
C.)  with  the  initial  mixing  conditions  such  as  the  velocity  ratio  or 
the  mass  flux  ratio.   The  results  will  be  discussed  in  Chapter  VI. 

Remarks  on  the  Turbulent  Transport  of  Mass  and  Energy 

In  early  analysis  of  mixing  problems,  it  was  often  assumed 
that  the  turbulent  Schmidt  and  Prandtl  numbers  were  unity;  an 
assumption  which  simplifies  the  governing  equations  considerably. 
However,  recent  experiments  indicate  that  the  Schmidt  number  may  differ 
significantly  from  unity.   Furthermore,  the  experimental  data  of  Forstall 
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and  Shapiro  [5]  show  that  the  Schmidt  number  remains  constant  at 

approximately  0.7  throughout  the  mixing  region,  so  that  the  eddy 

diffusivity  c,    is  merely  a  constant  times  the  eddy  viscosity  e    . 

For  gaseous  components  in  binary  mixing,  the  values  of  Sc  most 

frequently  cited  vary  betv/een  0.5  and  1.2.   In  the  present  investigation. 

the  turbulent  Schmidt  and  Prandtl  numbers  are  considered  as  parameters 

and  are  retained  constant  in  the  entire  mixing  region.   Using  suitable 

values  for  Sc  and  Pr  yields  the  values  for  e  ,  and  e,  from  calculated 

values  of  z    . 
m 

Hence,  with  the  choice  of  the  transport  mechanism,  the 
formulation  of  the  problem  is  completed.   The  equations  are  solved 
numerically  on  an  IBM  370/165  and  the  results  are  discussed  in 
Chapter  VI.   A  discussion  and  printout  of  the  computer  program  is 
presented  in  the  Appendix. 


CHAPTER  V 
EXPERIMENTAL  INVESTIGATION 

General  Considerations 

The  main  objective  of  the  experimental  investigation  was  to  examine 

the  nixing  of  parallel  two-dimensional  cof loving   supersonic  gas  streams. 

The  primary  emphasis  was  placed  on  obtaining  pressure,  velocity  and 

species  concentration  profiles  inside  the  mixing  region. 

The  secondary  stream  design  Mach  number  was  chosen  to  be  1.3 

for  the  following  reasons: 

i.   to  minimize  disturbances  at  the  nozzle  exit,  the  pressure  in 
the  two  streams  was  matched  by  presetting  settling  chamber 
conditions  in  each.   To  avoid  pressure  communication  back 
to  the  repective  chambers,  both  streams  had  to  be  supersonic, 
i.e.,  M  >  1.  Waves  due  to  supersonic  flow  at  Mach  1.3 
(if  any  existed)  would  tend  to  be  weak, 
ii.   an  average  run  time  of  approximately  thirty  seconds  was 

necessary  to  obtain  various  measurements.   Due  to  the  large 
number  of  runs  needed  to  accurately  define  the  specie 
profiles,  the  total  pressure,  in  the  secondary  settling  chamber 
would  have  to  be  as  low  as  possible  to  conserve  the  con- 
sumption of  commercial  bottled  gas.   High  Mach  numbers  in 
the  secondary  stream  would  have  forced  the  use  of  high 
total  pressures  to  match  exit  conditions.   Since  flow  rate 
is  directly  proportional  to  the  total  pressure,  a  design 
Mach  number  of  1.3  was  chosen  to  keep  the  secondary  stream 
mass  flow  rate  relatively  low. 
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iii.   the  contour  of  the  supersonic  portion  on  a  Mach  1.3  nozzle 

was  explicit  enough  to  ensure  relatively  error  free  machining 
The  primary  stream  design  Mach  numbers  of  1.3  and  2.0  were 
chosen  to  obtain  various  velocity  and  mass  flux  ratios  between  the  two 
streams.   A  summary  of  runs  with  various  configurations  and  test  con- 
ditions is  shown  in  Table  II. 


Experimental  Apparatus 

The  gas  dynamics  facilities  in  the  mechanical  engineering  de- 
partment were  modified  to  utilize  different  species  of  gas  for  the  mixing 
analysis  (Figure  4) .   A  block  diagram  of  the  system  is  also  presented 
in  Figure  5. 

The  existing  facility  consists  of  a  two-stage  positive  displace- 
ment type  compressor,  feeding  a  series  of  high  pressure  air  storage 
tanks.   Air  from  the  tanks  is  brought  to  the  laboratory  in  two  separate 
lines  to  a  pair  of  on-off  valves.   The  line  pressures  are  stepped  down 
to  the  desired  levels  by  a  pressure  regulator  in  each  stream  before  the 
air  enters  the  settling  chambers.   Inside  the  two  separate  settling 
chambers  are  a  series  of  flow  straighteners  and  dampening  screens.   Thus, 
fairly  uniform  streams  are  introduced  through  converging  sections  into 
the  test  section.   The  test  section  is  followed  by  a  diffuser  section; 
then  the  air  is  passed  through  a  sound  attenuator  and  exhausts  into  the 
atmosphere. 

This  system  was  modified  by  the  addition  of  a  storage  tank  that 
supplied  the  secondary  stream  with  either  argon  or  helium.   The 
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existing  dry  air  supply  system  was  utilized  for  the  primary  stream.   A 
two-dimensional  test  section  was  designed  that  permitted  the  use  of 
interchangeable  nozzle  blocks. 

Test  Section 


A  photograph  of  the  test  section  with  the  splitter  plate  and  a 
set  of  nozzle  blocks  installed  is  shown  in  Figure  6.   The  designed  test 
section  consisted  of  a  section  which  reduced  the  existing  dimensions  of 
the  system  to  the  desired  dimensions  of  8.5  by  0.5  inches  for  the  primary 
and  4.5  by  0.5  inches  for  the  secondary  stream.   This  section  was 
constructed  out  of  steel  and  had  a  secondary  function  of  supporting  most 
of  the  length  of  the  splitter  plate.   Erected  between  the  reducer  and 
the  diffuser  sections  was  the  main  frame  also  made  from  steel.   The 
upper  and  lower  portions  of  the  frame  were  used  to  position  and  secure 
the  primary  and  secondary  stream  nozzle  blocks.   Hard  neoprene  gasket 
material  was  used  to  seal  the  flanges.   The  frame  was  "sandwiched" 
between  two  one-inch  aluminum  side  plates  which  when  bolted  together  would 
give  a  test  section  width  of  one-half  inch.    Although  a  larger  test 
section  might  have  been  desirable,  the  flow  area  (thus  the  flow  rate) 
of  the  secondary  stream  was  the  governing  factor  in  the  test  section 
dimensions.   The  frame  and  the  side  plates  were  sealed  from  the  environ- 
ment by  linear  "0-rings".   The  side  plates  had  sections  cut  out  to  accom- 
odate optical  windows  and  the  static  pressure  plate. 

Since  having  the  two  streams  an  inf initesimally  small  distance 
apart  when  they  came  into  contact  with  each  other  was  physically  impos- 
sible, the  splitter  plate  was  machined  down  to  0.015  inches  at  the  tip. 
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A  thinner  plate  would  have  caused  strength  problems  since  the  splitter 
plate  had  to  support  the  force  due  to  the  pressure  difference  between 
the  streams. 

To  support  the  protruding  portion  of  the  splitter  plate,  a  1.25 
inch  wide  by  o.25  inch  deep  groove  was  machined  in  the  side  plates.  With 
the  splitter  plate  in  place,  epoxy  resin  was  poured  into  the  groove  and 
allowed  to  harden.   Since  the  splitter  plate  was  coated  with  silicone 
grease  before  this  operation,  it  was  easily  removed  when  the  resin  had 
hardened.   The  excess  resin  was  then  sanded  down  smooth  with  the  side 
plate  surface.   This  method  permitted  the  support  of  the  "odd-shaped" 
portion  of  the  splitter  plate.   The  snug  fit  also  served  as  a  seal 
between  the  primary  and  secondary  streams. 

Nozzle  Blocks 

The  contour  of  the  nozzle  blocks  guiding  the  subsonic  flow  was 
an  arbitrary  shape  which  permitted  smooth  transition  to  sonic  condi- 
tions at  the  minimum  area.   The  contour  providing  the  supersonic  flow 
was  determined  from  the  two-dimensional  method  of  characteristics. 
To  get  the  shortest  possible  test  section,  a  sharp-edged  throat  with 
a  single  wave  reflection  design  was  used.   No  allowance  for  the  boundary 
layer  was  made  in  the  nozzle  design;  however,  the  splitter  plate  had  a 
taper  of  0.007  inches  per  inch  at  the  straight  section.   This  made  up  for 
some  of  the  boundary  layer  accumulation  which  is  small  in  accelerated 
flows . 

The  nozzle  blocks  were  cut  out  of  0.5  inch  aluminum  plates, 
and  machined  to  the  desired  contour.   The  final  polishing  was  done  by 
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hand  using  fine  grained  emery  cloth.   Cushioned  tape  instead  of  O-rings 
was  placed  between  the  side  plates  and  the  nozzle  blocks  to  seal  the 
system. 

The  primary  stream  nozzles  had  a  throat  half-height  of  2  inches. 
Thus,  except  for  the  reflected  waves,  the  effect  of  the  wall  bounding 
the  primary  stream  could  be  neglected.   The  secondary  stream  nozzle  had 
a  throat  half-height  of  0.5  inches.   For  this  case  the  wall  effects 
could  not  be  ignored  and  limited  the  collection  of  data  to  the  downstream 
location  where  the  nixing  region  and  the  wall  boundary  layer  interacted 
with  each  other.   It  should  be  recognized  that  the  secondary  stream 
dimensions  were  determined  by  the  maximum  feasible  flow  rate  of  the 
stream. 

It  was  observed  during  calibration  runs  that  the  primary  and 
secondary  stream  Mach  1.3  nozzles  gave  surprisingly  clean  (shock-free) 
flows.   The  primary  stream  Mach  2.0  nozzle,  however,  did  display  some 
wave  patterns,  but  pressure  measurements  showed  these  waves  to  be  weak. 

The  actual  Mach  numbers  of  the  two  streams  were  checked  by 
three  methods:   1)   ratio  of  the  settling  chamber  total  pressure  to 
probe  total  pressure,  2)   ratio  of  the  local  static  pressure  to  probe 
total  pressure,  and  3)   measuring  the  wave  angles  on  the  schlieren 
photographs.   With  the  air  Mach  2.0  and  1.3  nozzles,  the  result  was  an 
average  Mach  number  of  1.97  and  1.28  respectively.   The  secondary  stream 
nozzle  (designed  for  M=1.3)  yielded  an  average  Mach  number  of  1.27. 
A  photograph  of  the  air  Mach  2.0  nozzle  block  is  presented  in  Figure  7. 
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Gas  Supply  and  Control  System 

The  primary  stream  utilized  air  from  the  existing  air  storage 
system  of  the  gas  dynamics  laboratory.   This  system  consists  of  a 
Worthington  two-stage  compressor  feeding  a  series  of  tanks  capable  of 
holding  approximately  420  cubic  feet  of  air  at  300  psi.   The  air  is 
passed  through  several  oil  and  water  traps  and  a  regenerative  type  gas 
dryer  before  it  goes  into  the  storage  tanks. 

The  secondary  stream  utilized  argon  (or  helium)  from  a  separate 
tank  capable  of  holding  approximately  50  cubic  feet  of  gas  at  150  psi 
pressure.   This  tank  itself  was  supplied  argon  (or  helium)  from  a  series  c 
commercial  bottled  gas  manifolded  together.   A  regulator  was  necessary 
to  reduce  the  commercial  pressure  from  a  maximum  of  approximately  2500 
psi  to  the  desired  150  psi. 

The  pressures  in  the  two  settling  chambers  were  controlled  by 
Fisher  pressure  regulators.   These  regulators  were  activated  by  remote 
control  with  a  single  switch,  and  if  necessary,  the  sequence  of  opera- 
tion could  be  staggered  with  the  bleed  valve  on  each  regulator. 
The  Schlieren  System 

A  schlieren  system  was  available  to  observe  the  behavior  of  the. 
flow  field.   A  xenon  lamp  was  used  as  the  monochromatic  light  source 
which  converged  on  a  16  inch  diameter  parabolic  mirror  through  a 
condensing  lens  and  a  knife  edge.   The  parallel  beam  of  light  reflected 
from  the  trirror  passed  through  the  test  section  and  was  reflected  off  of 
another  mirror  and  focused  on  a  knife  edge.   The  image  was  projected  on  a 
ground  glass  plate  and  by  the  use  of  a  Graflex  camera  polaroid  pictures 
of  this  image  were  taken.   Figures  8  and  9  show  two  typical  results 
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obtained.   The  shock  free  nature  of  the  flow  in  the  case  where  both 
streams  are  at  Mach  1.3  should  be  noted. 

Probe  Drive  Mechanism 

Due  to  the  small  dimensions  of  the  test  section,  it  was  necessary 
to  have  the  total  pressure  probe  very  close  to  the  primary  stream  wall 
during  the  start  of  the  run.   After  steady  flow  conditions  were  es- 
tablished in  the  test  section,  the  probe  had  to  be  passed  through  the 
free  stream  quickly  in  order  to  conserve  run  time.   It  then  had  to  be 
slowed  down  in  the  neighborhood  of  the  mixing  region  and  traverse  the 
mixing  region  at  a  relatively  slow  rate  to  obtain  an  accurate  pressure 
profile.   This  was  accomplished  by  the  use  of  a  variable  speed  reversible 
electric  motor  to  drive  the  probe. 

Eight  holes  were  drilled  at  various  intervals  on  top  of  the 
test  section  frame  and  the  nozzle  block  to  accomodate  the  probe  shaft. 
The  probe  shaft  was  attached  to  a  threaded  rod  which  in  turn  was  rotated 
by  the  electric  motor.   This  mechanism  was  mounted  on  top  of  the  test 
section  and  controlled  from  the  instrument  table.   Provisions  were  made 
to  accomodate  a  linear  potentiometer  and  a  pressure  transducer  en  the 
mechanism  frame.   Limit  switches  were  installed  to  automatically  stop  the 
probe  drive  once  it  reached  either  the  upper  or  lower  wall.   The  probe 
holes  not  in  use  were  closed  with  a  threaded  brass  plug  which  was 
screwed  in  until  the  tip  was  flush  with  the  nozzle  block  surface. 
Photographs  of  the  probe  drive  mechanism  and  the  mounting  on  the  test 
section  are  shown  in  Figures  10  and  11. 


yj 


Static  Pressure  Plate 

It  was  undesirable  to  obtain  static  pressure  measurements  with 
probes  in  such  a  small  test  section  since  the  presence  of  the  probe 
would  influence  the  results.   This  could  have  been  accomplished  by 
installing  pressure  taps  on  the  side  walls  if  it  were  not  for  the  schlieren 
windows.   An  alternative  was  to  machine  an  aluminum  plate,  the  exact 
size  of  the  windows,  which  would  be  interchangeable  with  one  of  the 
windows.   The  pressure  taps  could  then  be  drilled  in  this  plate.   This 
method  was  chosen  since  it  allowed  the  determination  of  static  pressure 
with  the  least  amount  of  external  disturbance  in  the  flow  field. 

A  total  of  ninety  holes  were  drilled  on  the  flow  field  side  with 
each  hole  having  a  diameter  of  0.030  inches.   These  holes  were  distributed 
among  13  axial  stations  ranging  from  0.25  inches  downstream  of  the  nozzle 
exit  plane  to  approximately  15  inches  downstream.   On  the  outside  of  the 
plate,  these  holes  were  enlarged  to  a  diameter  of  0.080  inches  so  that 
short  pieces  of  stainless  steel  tubing  could  be  pressed  in.   These  tubes 
were  sealed  at  the  bases  with  epoxy  resin.   Vinyl  tubing  was  used  to 
connect  the  stainless  steel  tubing  in  the  taps  to  the  monometer  board.   A 
photograph  of  the  static  pressure  plate  is  shown  in  Figure  12. 

Gas  Analysis  System 

Since  one  of  the  main  objectives  of  the  experimental  program 
was  to  obtain  species  concentration  profiles,  gas  samples  were  withdrawn 
from  the  flow  field  and  collected  in  a  series  of  vacuum  bottles.   These 
samples  were  subsequently  analyzed  on  a  Victoreen  Model  4000  Gas  chroma- 
tograph.   A  photograph  of  the  gas  chromatograph  together  with  the  columns 
used  is  shown  in  Figure  13. 
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The  fundamentals  of  gas  chromatography  can  be  found  in  reference 
[75].   The  chromatograph  measures  the  volumetric  concentration  of  each 
constituent  of  the  sample.   The  components  are  separated  when  passed 
through  a  column  consisting  of  a  length  of  stainless  steel  or  copper 
tubing  packed  with  a  solid  phase  such  as  charcoal.   Since  each  component 
progresses  through  the  column  at  different  rates,  the  travel  time 
(or  elution  time)  identifies  each  component  qualitatively.   Thermal 
conductivity  detectors  measure  the  quantity  of  each  of  the  separated 
gases  relative  to  the  carrier  gas  and  concentrations  are  printed  out 
on  a  strip  chart  recorder. 

Measuring  Devices 

Total  Pressure  Measurements 

Total  pressure  profiles  in  the  mixing  region  were  obtained  by 
introducing  a  probe  into  the  stream.   The  pressure  registered  via  the 
probe  was  transformed  into  an  electrical  output  using  a  MB  Electronics 
Model  151-BAA-l  pressure  transducer.   Power  was  supplied  to  the  trans- 
ducer by  a  CEC-3-140  DC  power  supply.   The  output  was  recorded  on  a 
CEC-5124A  20  channel  recording  oscillograph.   This  strip  chart  recorder 
made  traces  on  light-sensitive  tape  which  was  6  inches  wide.   Thus,  the 
voltage  output  from  the  transducer  had  to  be  scaled  down  by  means  of 
an  external  attenuation  circuit  so  that  a  full  scale  deflection  registered 
180  psia  pressure.   With  this  calibration,  probe  pressures  could  be 
read  to  within  +0.75  psia. 

The  probes  had  to  be  small  enough  so  that  the  least  amount  of 
disturbance  would  be  introduced  into  the  flow  field  and  thus  give 
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accurate  pressure  readings.   Yet  they  had  to  be  strong  enough  to  withstand 
the  bending  moments  due  to  high  speed  flow.   For  this  purpose  stain- 
less steel  tubing  of  0.060  inches  OD  (0.036  inch  ID)  was  used.   A 
length  of  this  tubing  was  bent  at  a  90  degree  angle  and  welded  into  a 
short  piece  of  0.25  inch  diameter  stainless  tubing  which  in  turn  screwed 
into  the  prcbe  shaft  on  the  drive  mechanism.   The  length  of  the 
"sting"  was  determined  by  the  location  of  the  probe  holes  in  the  test 
section  frame  relative  to  the  static  pressure  taps  since  total  pressure 
profiles  were  needed  at  the  point  where  static  pressure  data  were  taken. 
It  was  found  that  three  different  "sting"  lengths  of  0.2,  0.5  and  0.75 
inches  were  needed.   After  the  desired  lengths  were  cut,  a  0.005  inch 
thick  shim  stock  was  inserted  into  the  end  of  the  sting  and  compressed 
to  form  a  slit  0.005  by  0.040  inches.   This  design  yielded  satisfactory 
results  in  terms  of  accuracy  and  strength.   A  photograph  of  the  probes 
is  shown  in  Figure  14. 

The  pressures  of  the  two  settling  chambers  were  monitored  both 
visually  on  pressure  gauges  and  also  on  the  recording  tape  by  the  use 
of  two  Giannini  Model  46139  pressure  transducers. 

Total  Temperature  Measurements 

Total  temperatures  were  monitored  only  in  the  two  settling  cham- 
bers using  chrome 1-alumel  alloy  thermocouples,  the  output  of  which  was 
recorded  on  the  recording  tape.   The  reference  junction  was  held  at 
32°F  by  immersing  it  in  an  ice  bath.   A  full  scale  deflection  of  six 
inches  on  the  recording  tape  corresponded  to  temperature  readings 
between  32-92  °F,  which  was  the  temperature  range  of  interest. 
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Static  Pressure  Measurements 

A  thirty-tube  illuminated  mercury  manometer  was  used  to  measure 
static  pressures.   A  pressure  differential  of  approximately  60  inches  of 
mercury  could  be  measured  on  this  manometer.   By  using  atmospheric  pressure 
as  reference  and  setting  the  zero  at  the  mid-point,  pressure  measure- 
ments in  the  range  of  0-30  psia  were  obtained. 

Since  there  were  ninety  pressure  taps  available  and  only  thirty 
manometer  tubes,  a  manifolding  system  had  to  be  devised.   This  was  accom- 
plished by  connecting  each  manometer  tube  to  a  six  inch  length  of  brass 
pipe  at  one  end  and  welding  shut  the  other.   Three  brass  stop-cock  valves 
were  mounted  in  each  piece  of  pipe  and  a  pressure  tap  connected  to  each 
valve.   Thirty  of  these  manifolds  were  mounted  on  the  manometer  board. 
Vinyl  tubing  was  used  to  make  the  connections  and  each  connection  was 
sealed  with  enamel  paint  (Glyptal).   No  leakage  problems  were  encountered. 
Thus,  each  manometer  tube  was  capable  of  reading  one  of  three  pressure 
taps  depending  on  which  valve  was  turned  on. 

A  Graflex  camera  using  4  by  5  Polaroid  plate  film  was  used 
to  record  the  pressure  measurements.   From  the  photographs,  the  pressures 
could  be  read  to  within  0.1  inches  of  mercury  allowing  pressure  measure- 
ments to  within  +0*05  psia. 

Concentration  Measurements 


1)  Probe  design 

In  order  to  reduce  the  number  of  runs  necessary  for  the  accurate 
determination  of  concentration  profiles,  gas  sample  rakes  made  up  of 
three  probes  each  were  designed.   The  tips  of  the  probes  were  made  from 
0.040  inch  OD  by  0.009  inch  wall  stainless  steel  tubing.   The  tips  were 
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then  immediately  expanded  to  0.040  inch  ID  to  prevent  the  flow  inside 
the  probe  from  choking.   As  was  done  in   the  case  of  the  pressure  probes, 
the  tubing  was  bent  at  a  90  degree  angle  and  mounted  inside  a  0.25 
inch  diameter  by  1.0  inch  long  stainless  steel  tubing  which  in  turn 
screwed  into  the  probe  shaft  on  the  drive  mechanism.   The  nominal 
distance  between  the  probe  tips  was  0.1  inches.   Each  of  the  three  probes 
was  connected  to  a  2.5  cubic  inch  volume  evacuated  sample  bottles  by 
means  of  vinyl  tubing.   A  photograph  of  the  gas  sampling  rakes  and  a 
sample  collection  bottle  is  shown  in  Figure  15. 
2)  Chrcmatograms 

The  column  used  in  the  gas  chromatograph  was  a  6-foot  length 
of  Varian  5-A  molecular  sieve.   The  column  was  conditioned  by  drying 
it  in  the  chromatograph  oven  set  at  750  CF  for  24  hours.   During  the 
conditioning  time  the  carrier  gas  was  allowed  to  flow  through  the  column 
at  a  rate  of  40  ml/minute.  A  Varian  Model  02-001126-00  thermal  con- 
ductivity cell  using  two  pairs  of  30  ohm  tungsten- rhenium  filaments  was 
used  to  measure  the  amount  of  each  constituent  in  the  mixture.   The 
output  from  the  conductivity  cell  was  recorded  by  a  Honeywell  Model 
Electronik-194  strip  chart  recorder. 

In  the  Series  I  tests  where  argon-air  mixtures  were  being  analyzed, 
helium  was  used  as  the  carrier  gas ;  argon  was  used  as  the  carrier  gas  to 
analyze  helium-air  mixtures.   For  both  series  of  tests,  the  injection 
ports,  the  column  and  the  thermal  conductivity  cell  were  maintained  at  86  °F. 


1The  spacing  between  the  sample  probes  was  determined  from  the 
recommendations  of  reference  [17].   With  the  above  design  no  inter- 
ference problems  were  encountered. 
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The  Series  II  tests  of  helium-air  mixtures  presented  no  problems 
since  the  molecular  sieve  column  separated  oxygen,  nitrogen  and  helium  into 
distinct  bands  and  the  mass  fractions  could  be  calculated  from  peak 
areas.  However,  argon  cannot  be  separated  from  oxygen  when  column 
temperatures  are  above  approximately  -95  °F  [76,  77],   If  an  acetone- 
dry  ice  bath  was  used  to  attain  this  temperature  then  nitrogen  would  not 
be  eluted  and  thus  be  irreversibly  adsorbed  in  the  column.   Another 
alternative  was  to  separate  nitrogen  from  the  oxygen-argon  mixture  at 
room  temperature,  then  immerse  the  column  in  the  dry  ice-acetone  bath 
and  inject  the  sample  again.   For  a  single  chromatogram  of  a  mixture,  the 
turn-around  time  using  this  technique  was  estimated  to  be  over  30  minutes. 
It  was  necessary  to  have  three  or  more  chromatograms  of  the  same  mixture 
to  obtain  a  statistical  average  of  the  mass  fractions  of  the  constituents. 
This,  together  with  the  large  number  of  samples  collected,  made  the  above 
method  impractical. 

An  indirect  method  of  calibration  was  devised  to  avoid  this  problem. 
It  was  assumed  that  air  behaved  as  a  single  specie^.   As  a  measure  on 
the  validity  of  this  assumption,  self-diffusion  coefficients  obtained  by 
kinetic  theory  considerations  [78]  were  examined.   It  was  found  that  in 
the  temperature  and  pressure  range  of  interest,  the  self-diffusion 


Air  is  conventionally  treated  as  a  single  component  in  evalu- 
ation of  transport  properties  for  low  density  systems  (i.e.  ,  pressures 
on  the  order  of  one  atmosphere)  [62].   Example  calculations  of  binary 
and  ternary  diffusion  in  air  (i.e.  ,  considering  it  as  a  single  specie 
and  as  a  mixture)  may  be  found  in  reference  [62].   As  is  usually  expected, 
the  two  methods  are  in  good  agreement. 
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coefficients  were  within  5  percent  of  each  other.   Thus,  for  diffusion 
purposes,  the  nitrogen  and  oxygen  molecules  were  practically  indistingu- 
ishable.  Then,  if  the  mass  fraction  of  nitrogen  in  the  mixture  was  known, 
the  mass  fraction  of  oxygen  could  be  computed  since  for  every  0.79 
moles  of  nitrogen  there  are  0.21  moles  of  oxygen.   Knowing  the  mass 
fractions  of  two  of  the  three  constituents,  the  third  could  be  deduced. 

Exact  amounts  of  pure  argon,  oxygen  and  nitrogen  were  injected 
into  the  column.   The  values  of  the  peak  areas,  calculated  by  the  method 
of  triangulation,  were  then  plotted  against  the  known  weight  injected. 
As  expected  the  curves  were  linear  and  passed  through  the  origin  (no 
sample,  no  response). 

When  an  exact  amount  of  unknown  sample  was  injected  into  the 
column  at  room  temperature,  two  peaks  would  appear  on  the  chromatogram; 
one  of  pure  nitrogen  and  one  of  the  argon-oxygen  mixture.   Using  the 
calibration  curves  the  weight  of  nitrogen  could  be  determined  from  its 
peak  area.   Then  a  simple  ratio  would  yield  the  weight  of  oxygen  present 
in  the  mixture.   The  calibration  curves  would  again  be  used,  this  time 
somewhat  in  reverse,  to  obtain  the  peak  area  corresponding  to  this  weight. 
This  area  would  then  be  subtracted  from  the  "compound"  argon-oxygen  peak 
area  to  get  the  argon  peak  area  and  therefore  its  weight.   Knowing  the 
weights  of  each  of  the  constituents  of  the  mixture,  the  mass  fractions 
could  be  calculated. 

It  was  realized  from  the  start  that  errors  would  be  magnified 
when  the  concentration  of  argon  in  the  mixture  fell  below  approximately 
10  percent,  i.e.,  samples  collected  from  cross-stream  locations  close 
to  the  primary  stream.   The  accuracy  of  the  overall  data,  discussed  in 


66 


the  next  section,  showed  this  error  to  be  at  an  acceptable  level. 

To  keep  all  other  sources  of  error  at  a  minimum,  the  following 
precautions  were  taken: 

i.   calibrations  were  frequently  checked. 
ii.   injection  port,  column  and  detector  cell  temperatures 

together  with  the  carrier  gas  flow  rate  were  maintained 
at  the  same  levels  as  used  in  calibration  runs, 
iii.   fast  recorder  chart  speed  was  used  to  make  peak  width 
measurements  more  accurate, 
iv.   signal  attenuations  were  adjusted  to  obtain  full  scale  peaks 
so  that  peak  height  measurements  were  accurate, 
v.   the  mximum  allowable  filament  current  was  used  in  the  detector 
cell  to  increase  the  overall  accuracy  of  the  chromatograms. 
vi.   number  of  chromatograms  per  mixture  were  increased  to  get 
better  statistical  values. 

Testing  Procedure 

Since  all  of  the  measurements,  namely  static  and  total  pressure, 
concentration  and  schlieren  photographs,  could  not  be  made  during  the 
course  of  one  run,  it  was  necessary  to  make  different  sets  of  runs.   The 
data  collection  technique,  repeated  for  each  set  of  nozzle  blocks,  is 
summarized  below: 

1)  The  settling  chamber  total  pressures  necessary  to  match 
static  pressures  in  the  exit,  plane  of  the  nozzles  were  calculated. 
Using  air  in  both  streams  the  regulators  were  adjusted  to  yield  these 
pressures.   The  secondary  stream  air  line  was  then  shut  off  and  the  static 
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pressure  plate  was  installed  in  the  test  section.   With  the  secondary 
stream  now  utilizing  argon  (or  helium) ,  further  adjustments  were  made 
in  the  chamber  total  pressures  by  observing  the  exit  plane  static  pressures, 
The  static  pressure  plate  was  removed  and  the  schlieren  window  was  in- 
stalled, after  which  the  flow  field  was  observed  for  "cleanliness", 
i.e.,  the  presence  of  unwanted  expansion  or  compression  waves.   The 
repeatability  of  chamber  pressures  from  run  to  run  was  within  3.5 
percent  for  the  primary  and  1.5  percent  for  the  secondary.   Better 
repeatability  (i.e.,  1  percent)  was  obtained  with  the  primary  stream 
when  Ma eh  1.3  nozzle  block  was  used. 

2)  After  a  satisfactory  flow  field  was  established,  a  schlieren 
photograph  was  made.   The  static  pressure  plate  was  replaced  and  three 
runs  made  to  determine  the  static  pressure  distribution.   In  between  the 
runs  the  appropriate  manifold  valves  were  turned  on. 

3)  The  total  pressure  runs  were  made  with  the  windows  back  in 
place  so  that  the  shock  pattern  due  to  the  presence  of  the  probe  in  the 
flow  field  could  be  viewed  on  the  schlieren  screen.   The  number  of 

runs  at  each  axial  location  depended  on  repeatability  and  the  quality  of 
the  traces.   As  the  pressure  probe  was  traversing  the  flow  field,  the 
exact  location  of  the  probe  tip  had  to  be  known.   To  accomplish  this  a 
linear  potentiometer  was  connected  to  the  travelling  probe  mechanism. 
Before  each  run,  position  calibration  was  done  by  getting  a  trace  on 
the  recorder  tape  while  the  probe  tip  was  located  at  the  lower  wall. 
Thus,  knowing  a  reference  position  and  measuring  the  displacement  of  the 
potentiometer  output  trace  obtained  during  a  run,  the  position  of  the 
probe  and  the  corresponding  pressure  at  that  point  could  be  obtained. 
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At  the  begining  of  each  run  the  probe  was  located  next  to  the  primary 
stream  wall.   The  "blow-down"  was  started  and  allowed  to  reach  a  steady 
state  after  which  the  probe  mechanism  was  activated.   The  probe  was  allowed 
to  "sweep"  the  entire  flow  field,  but  as  it  approached  the  lower  wall, 
the  direction  was  reversed  and  the  speed  reduced  so  as  to  obtain  a 
"fine- trace"  through  the  mixing  region. 

4)  A  vacuum  pump  was  connected  to  a  manifold  with  three  outlets; 
one  outlet  was  connected  to  a  30  inch  vertical  mercury  manometer,  the 
second  to  the  injection  port  of  the  gas  chromatograph  and  the  third  to 
the  sample  bottle.   Before  each  sample  collection  run,  the  sample  bottles 
were  evacuated,  after  which  the  system  was  purged  with  carrier  gas  so  as 
to  minimize  the  concentration  of  any  possible  residual  sample  from  the 
previous  run.   Position  calibration  was  done  in  the  same  manner  as  was 
done  in  total  pressure  measurements,  by  obtaining  a  trace  with  the 
outermost  probe  on  the  lower  wall.   The  probe  rake  was  then  moved  to  the 
location  where  sample  collection  was  desired  by  activating  the  drive 
mechanism.   This  location  was  determined  from  schlieren  photographs; 
that  is,  most  of  the  samples  were  collected  from  inside  the  mixing 
region  with  only  a  couple  of  measurements  in  the  free  streams.   The 
"blow-down"  was  started  and  the  flow  of  sample  was  established  through  the 
probes.   After  steady  flow  conditions  prevailed,  the  probes  were  connected 
to  the  sample  bottles  and  the  bottle  valves  were  opened.   It  was  observed 
that  spproximately  20  seconds  of  run  time  was  required  to  obtain  an 
"adequate1"  quantity  of  sample.   After  each  run  the  samples  were  normalized 


1 
"adequate"  quantity  was  determined  by  trial  and  error  to  be  a 
sample  at  approximately  1/3  to  1/2  atmosphere  pressure.   This  yielded  a 
sample  of  high  enough  concentration  after  normalization  to  one  atmosphere 
pressure  for  the  required  number  of  chromatograms  (i.e.,  3-6). 
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to  one  atmosphere  pressure  by  the  addition  of  carrier  gas  into  the  bottles. 
Then  1  ml.  of  each  sample  was  injected  into  the  gas  chromatograph  by 
means  of  a  micro-volume  gas  sampling  valve  which  is  an  integral  part  of  the 
gas  chromatograph.   A  slight  vacuum  had  to  be  applied  to  the  exhaust  port 
of  the  sampling  valve  to  "suck"  in  the  sample  which  was  at  one  atmosphere 
pressure. 

Data  Reduction 


Two  methods  of  reducing  the  data  were  considered.   The  first 
method  was  to  assume  uniform  static  pressure  equal  to  an  average  test 
section  pressure  throughout  the  flow  field.   The  second  method  was  to  use 
the  average  static  pressure  at  each  downstream  station  after  making  sure 
that  the  cross-stream  variation  was  less  than  +  5  percent.   Either 
method  could  be  used  to  reduce  the  data  of  Series  I-B  and  II-B  tests 
with  practically  the  same  results  since  the  streamwise  static  pressure 
variation  is  very  small  (Figures  17  and  19).   On  the  other  hand,  as  may 
be  seen  in  Figures  16  and  18,  some  static  pressure  variations  were 
observed  in  the  Series  I-A  and  II-A  tests.   Both  methods  were  used  to 
reduce  the  data  of  Series  I-A  and  II-A  tests.   It  was  observed  that  the 
two  methods  yielded  reduced  velocity  data  within  2-3  percent  of  each 
other,  provided  the  streamwise  static  pressure  variations  were  within 
6-8  percent  of  the  exit  plane  value.   Hence,  the  results  of  the  second 
method  of  data  reduction  are  presented  because  retention  of  a  true 
representation  of  conditions  in  the  mixing  region  is  desirable.   It  is 
noted,  however,  that  an  average  test  section  static  pressure  was  used  in 
the  mixing  analysis  of  Chapter  III. 
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Since  the  two  streams  had  different  specific  heat  ratios,  con- 
centration values  were  needed  in  conjunction  with  the  pressure  measure- 
ments.  The  total  pressure  data  obtained  re.rein  the  form  of  a  smooth 
trace  from  the  recorder,  whereas  concentration  measurements  were  points 
spread  throughout  the  mixing  region.   Total  pressure  and  concentration 
measurements  corresponding  to  a  specific  cross-stream  location  were 
obtained  by  plotting  the  total  pressure  and  concentration  profiles  on 
the  same  graph,  then  drawing  a  smooth  curve  through  the  concentration 
values;  thus,  the  species  concentration  and  the  total  pressure  corres- 
ponding to  a  specific  cross-stream  location  were  available. 

The  following  equations  were  used  to  obtain  the  average 
specific  heat  ratios  of  the  mixtures: 

C  =  I   a.C  .  (5.1) 

P   J  ipi 

where   C  =  average  C  of  the  mixture  (btu/lb-R) 
P  P 

a.  ■  mass  fraction  of  species i  (lb  i/lb  mixture) 

C  .  =  C   of  speciesi  (Btu/lb-R) 
pi    P 

W-  I/I  ^  (5.2) 

i  Wi 

where    W  =  average  molecular  weight  of  mixture 

W.  =  molecular  weight  of  speciesi 

R=^  (5.3) 

W 

where   R  =  gas  constant  of  mixture  (Btu/lb-R) 

&  =  universal  gas  constant  (Btu/lb  mole-R) 

k  =  _  P  _  (5.4) 

C   -  R 

P 
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where   k  =  average  specific  heat  ratio  of  mixture. 

Velocities  were  computed  by  first  determining  the  Mach  numbers 
through  the  Rayleigh  pitot  formula  [79]: 


p       ui   9  k/k_1   o,    o   i  i  lk~1 
_5y.=  (k±iM2)      /(2k_  2_k-l 

P     k  2  "/      ,K   k+1      k+1  ;  *   ; 


where   P   =  the  measured  probe  pressure  at  (x,y) 

P  =  the  measured  average  static  pressure  at  (x) 
k  =  the  average  specific  heat  ratio  a  (x,y) 
M  =  the  local  Mach  number. 
Since  local  values  of  static  and  total  pressures  in  conjunction  with 
mass  fractions  are  used,  the  utilization  of  Equation  (5.5)  for  mixtures 
is  justified.   In  other  words,  no  attempt  is  being  made  to  relate  any 
of  the  quantities  along  streamlines  to  the  undisturbed  portion  of  the 
flow. 

Although  the  distribution  of  total  temperature  across  the 
mixing  zone  is  non-uniform  even  for  streams  of  equal  total  temperatures, 
this  non-uniformity  is  small  if  the  Prandtl  number  is  close  to  unity. 
For  example,  in  air-air  mixing,  the  variation  of  the  total  temperature 
is  only  about  0.1  percent  for  a  Mach  number  of  unity  [1].   In  his 
investigation  of  the  supersonic  mixing  of  hydrogen  and  air,  Morgenthaler 
[80]  observed  that  typical  experimental  profiles  at  Mach  2.0  indicated 
a  3  percent  variation  in  total  temperature. 

On  this  basis  the  total  temperature  was  assumed  to  be  constant 
through  the  mixing  region  and  the  average  of  the  total  temperatures  of 
the  two  streams  was  used  for  T  .   The  total  temperatures  of  the  two 
streams  were  never  more  than  10°  apart  with  the  average  value  being 
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approximately  535  R.   Thus,  static  temperature  profiles  were  determined 
through  the  relation 

T  -  T  /(l  +  -^M2)  (5.6) 

0  l 

For  a  mixture  of  n  species,  the  local  mass  average  velocity  U 

is  defined  as  [62] 

n 

1  P-u- 

1-1  X   1    ? 

U  =  — =   )   cx.U.  (5.7) 

n         '   i  i 

r  1=1 

It  is  noted  that  U  is  the  velocity  one  would  measure  by  means  cf  a 
pitot  tube  (i.e.,  incompressible  flow)  and  corresponds  to  the  velocity 
as  used  for  pure  fluids.   The  local  velocity  of  each  species  was 
calculated  using  the  adiabatic  flow  equation 


Ui  =  V2Cpi(To  "  T)  (5'8) 

Finally,  densities  were  computed  from  the  perfect  gas  law 


P  -  ~  (5.9) 

RT 

The  assumption  of  argon,  helium,  oxygen  and  nitrogen  being  perfect 

gases  was  valid  since  the  pressures  were  much  less  than  the  critical 

pressure  and  the  temperatures  much  greater  than  the  critical 

temperature  for  all  species  during  all  test  conditions. 

It  was  desired  to  relate  the  reduced  velocity  and  temperature 

data  to  a  set  of  average  initial  conditions  for  each  test  configuration 

(i.e.,  an  average  primary  stream  velocity  and  temperature  together 

with  an  average  secondary  stream  velocity  and  temperature).   This  was 
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necessary  because  the  theoretical  analysis  in  Chapter  III  required 
initial  profiles  of  velocity  and  temperature  and  with  these  profiles 
the  mixing  program  "marched"  downstream.   Since  the  experimental 
system  was  not  perfectly  repeatable,  each  set  of  initial  conditions 
varied  somewhat  from  run  to  run.   In  addition  to  the  above,  in  some 
of  the  experimental  configurations  as  discussed  earlier,  there  were 
slight  variations  in  the  static  pressure  which  affected  the  velocity 
profiles.   Eggers  and  Torrence  [81],  in  their  experimental  investigation 
of  compressible  air  jets  encountered  similar  problems.   They  compensated 
for  the  above  stated  variations  by  suggesting  a  velocity  modification 

of  the  type 

V  -  U     V  -  U 

§_  =  _£ §x_ 

U   -  U    U   -  U  \->-±vj 

p    s    px    sx 

where   V  =  the  new  local  velocity  modified  for  pressure  and  free 
stream  deviations 
U  =  the  average  primary  stream  velocity  at  the  nozzle  exit 

plane 
U  -  the  average  secondary  stream  velocity  at  the  nozzle  exit 

plane 
V  =  the  local  velocity  modified  for  static  pressure  changes 
U   =  the  primary  stream  velocity  modified  for  static  pressure 
changes 
the  secc 
changes 
The  above  modification  was  adopted  for  the  present  analysis  together 
with  a  temperature  modification  of  the  same  type 
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s_  _  p    sx 

p     s     px     SX 


where   T  =  the  new  local  temperature  (modified) 

T  =  the  average  primary  stream  temperature  at  the  nozzle 

exit  plane 
T  =  the  average  secondary  stream  temperature  at  the  nozzle 

exit  plane 

T'  =  the  local  temperature 
P 

T   =  the  local  primary  stream  temperature 

T   =  the  local  secondary  stream  temperature. 

It  should  be  recognized  that  these  modifications  were  adopted 
for  the  sake  of  consistency  between  the  experimental  data  and  the 
theoretically  predicted  profiles.   Since  both  the  experimental  and 
theoretical  results  are  presented  in  the  form  of  excess  velocity 
profiles  (discussed  in  the  next  chapter)  there  are  no  consequences  due 
to  these  modifications. 

Accuracy  of  Results 

Based  on  the  chart  and  photograph  resolutions  of  recorded  data, 
repeatability  cf  runs  and  calibration,  the  estimated  accuracy  of  the 
measurements  are: 

Static  pressures  1  0.1  psia 

Probe  pressures  t   0. 75  psia 

Total  temperatures  ±  10° 

Frobe  position  1  0 . 025  inches . 

The  test  section  static  pressures  are  on  the  order  of  10  psia; 
the  probe  pressures  range  from  35-70  psia;  the  temperature  range  of 
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interest  is  on  the  order  of  535  R  and  measurements  are  made  in  a  region 
of  approximately  2  inches.   Thus,  in  terms  of  percentage  errors: 

Static  pressures  *  1% 

Probe  pressures  -   2.5% 

Total  temperatures  i  2% 

Probe  position i  1.5% 

With  the  above  values  and  the  equations  used  to  reduce  the 
data,  it  is  estimated  that  the  velocity  data  are  accurate  to  within 
3  percent.   This  estimation  does  not  include  any  uncertainty  due  to 
concentration  measurements. 

One  means  of  assessing  the  overall  accuracy  of  the  data  is 
to  apply  the  principle  of  conservation  of  mass  to  the  secondary 
stream.   This  would  also  indicate  the  degree  of  accuracy  of  the 
concentration  measurements.   The  following  equation  was  evaluated 
numerically  for  each  axial  station  at  which  data  were  taken: 


mc  =  J  a.pv 


m_  =   a,.pudA  (5.12) 

A 
where   pu  =  local  mass  flow  per  unit  area  evaluated  from  experimental 

data 

a.  -  local  mass  fraction  of  the  secondary  stream  constituent 

A  =  area  of  the  flow  field  over  which  a.  is  nonzero. 

i 

If  the  data  were  correct,  m  would  be  equal  to  the  secondary 
stream  flow  rate  which  can  be  approximated  from  settling  chamber 
conditions  (for  uniform,  one-dimensional,  isentropic  flow) 

.    /r 

—  T±=   const.  (5.13) 

t    o 
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where  the  constant  depends  on  whether  argon  or  helium  is  being  used. 

Differences  between  rn  and  m  are  due  to  experimental  error.   A 
c 

comparison  between  these  two  mass  flow  quantities  is  an  essential 
criterion  in  assessing  the  accuracy  of  concentration  measurements 
because  of  uncertainties  in  obtaining  representative  samples  from 
flowing  streams. 

The  application  of  this  criterion  to  the  data  presented  herein 
is  reported  in  Table  III.   Large  errors  in  concentration  measurements 
taken  from  binary  streams  may  occur;  the  probe  design,  sampling  technique 
and  the  local  turbulence  level  in  the  flow  field  have  a  significant 
effect  upon  the  results  [35,81].   The  actual  physical  mechanism  which 
causes  unrepresentative  sample  collection  is  not  known,  but  satisfactory 
results  were  obtained  with  the  probe  and  sampling  technique  used  in 
this  investigation.   As  can  be  seen  in  Table  III,  the  overall  error 
is  less  than  10  percent  for  all  cases,  and  75  percent  of  the  cases 
have  an  error  of  6  percent  or  less.   Similar  sampling  problems  were 
found  in  references  [35 , 81, 82 ,83 ,84 ,85]  where  errors  of  up  to  ±  25 
percent  were  encountered.   Therefore,  it  is  concluded  that  the  accuracy 
of  the  data  is  well  within  acceptable  limits.   Furthermore,  it  is 
deemed  that  the  error  involved  in  the  method  of  determining  argon  mass 
fractions  in  argon-air  mixtures  is  negligible  as  compared  to  the 
uncertainties  of  the  samples  themselves. 


CHAPTER  VI 
DISCUSSION  OF  THEORETICAL  AND  EXPERIMENTAL  RESULTS 

Schlieren  Photographs 

Schlieren  photographs  were  made  to  observe  the  quality  of  the 
supersonic  flow  for  each  of  the  four  test  configurations.   Figure  8 
shows  the  flow  field  with  primary  stream  of  air  at  a  Mach  number  of 
approximately  2.0  and  the  secondary  stream  of  argon  at  a  Mach  number 
of  about  1.3.   Here,  the  wave  patterns  are  distinct  and  although  the 
mixing  zone  is  not  too  clearly  visible,  the  waves  may  be  observed 
to  bend  as  they  pass  through  the  mixing  zone.   The  "left-running"  wave 
emanating  from  the  left  middle  of  the  photograph  is  the  typical 
"lip  shock"  resulting  from  two  supersonic  streams  coming  into  contact 
with  each  other.   Static  and  total  pressure  measurements  confirmed 
these  waves  as  being  weak  and  the  local  wave  angles  in  the  primary 
stream  indicated  a  Mach  number  on  the  order  of  1.97. 

There  is  considerable  difference  between  the  flow  fields 
depicted  in  Figures  8  and  9.   In  Figure  9  both  streams  are  at  a  Mach 
number  of  1.3  again  with  air  in  the  primary  and  argon  in  the  secondary 
streams.   The  possibility  of  the  flow  being  subsonic  due  to  the 
absence  (or  non-visibility)  of  the  "lip  wave"  and  other  waves  was 
discarded  with  total  pressure  measurements.   When  a  total  pressure 
probe  was  injected  into  the  flow  field,  weak  oblique  shock  waves 
were  also  observed  around  the  probe  tip.   Due  to  the  cleanness  of  the 
flow,  the  mixing  zone  is  more  distinct  than  in  Figure  8. 

Similar  qualitative  results  were  obtained  when  helium  was 
utilized  in  the  secondary  stream.   The  small  width  of  the  test  section 
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was  a  factor  in  the  quality  of  the  schlieren  photographs  since  the 
quality  of  the  image  is  a  function  of  the  width  over  which  the 
initially  parallel  light  beams  are  diffracted. 

Static  Pressure  Variations 


The  mixing  analysis  of  this  study  involved  the  use  of  the 
boundary  layer  form  of  the  conservation  equations,  from  which  it  was 
deduced  that  the  transverse  pressure  gradient  (9P/Dy)  was  negligible. 
It  was  further  assumed  that  the  streamwise  pressure  gradient  (cP/3x) 
could  also  be  neglected.   The  validity  of  these  assumptions  are  now 
analyzed  in  light  of  the  experimental  data  obtained.   The  static 
pressure  measurements  are  plotted  in  Figures  16  through  27  with  the 
average  static  pressure  at  the  exit  plane  of  the  nozzles  used  as  a 
reference  pressure  and  with  the  physical  coordinates  non-dimensionalized 
with  respect  to  the  exit  height  of  the  secondary  stream,  i.e.,  slot 
height. 

For  all  four  test  conditions,  the  streamwise  pressure  distri- 
bution is  plotted  along  three  lateral  locations;  along  the  plane  of 
the  splitter  plate  (y/a  =  1.0),  and  one-half  slot  height  above  and 
below  the  plane  of  the  splitter  plate  (y/a  =  1.5  and  y/a  =  0.5).   It  may 
be  observed  in  Figures  16  through  19  that  the  static  pressure  increases 
monotonically  after  a  certain  axial  location.   This  is  a  typical 
characteristic  of  confined  (ducted)  flows.   The  axial  location  at  which 
this  steady  increase  is  observed  usually  corresponded  to  approximately 
the  downstream  location  where  the  mixing  zone  interacted  with  the  wall 
boundary  layer.   UTien  both  streams  are  at  Mach  1.3,  Series  I-B  and 
II-B  tests  (Figures  17  and  19) ,  the  static  pressure  variation  is  within 
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2-3  percent  of  the  average  exit  plane  value  for  about  12  slot  heights. 
In  Series  I-A  and  II-A  tests,  where  the  primary  stream  of  air  is  at 
Mach  2.0  and  the  secondary  stream  of  argon  or  helium  is  at  Mach  1.3 
(Figures  16  and  18) ,  higher  variations  of  up  to  10  percent  for  regions 
within  12  slot  heights  are  observed.   When  the  streamwise  static 
pressure  variations  and  the  schlieren  results  are  considered  together, 
it  is  concluded  that,  for  the  configurations  involving  the  Mach  2.0 
nozzle  block,  the  variation  of  the  static  pressure  can  be  attributed 
to  the  weak  waves  present  in  the  flow  field. 

A  further  check  on  the  magnitude  of  the  streamwise  pressure 

gradient  was  also  made  by  comparing  it  to  one  of  the  convective  terms 

?  1 

(i.e.,  pu")  in  the  momentum  equation.    The  results  are  presented 

below  in  terms  of  the  parameter  9 ,  where  6  is  defined  as 


|  »*2*y'h\    Pdy 


x/a  9 

1.0  55 

3.0  42 

5.0  39 

7.0  36 

9.0  32 

The  above  values  are  from  the  data  of  Series  II-A  tests.   Similar 

results  were  obtained  with  the  Series  I-A  test  data. 

Therefore,  if  only  the  region  upstream  of  the  point  of 


These  values  were  obtained  by  getting  intermediate  printed  out- 
put from  the  computer  program. 
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interaction  of  Che  mixing  region  and  the  confining  walls  is  considered, 
the  assumption  of  constant  axial  pressure  distribution  is  well  justified 
for  one  configuration  and  at  least  acceptable  for  the  other. 

The  transverse  pressure  variations  for  the  Series  I  and  II 
test  at  various  axial  locations  are  presented  in  Figures  20  through 
27.   For  Series  I-B  tests  involving  both  streams  at  the  same  Mach 
number  of  1.3,  the  variation  is  within  -  5  percent  for  all  stations  up 
to  x/a  =  12.   The  same  is  essentially  true  for  Series  I-A  tests  except 
that  at  about  x/a  =  12  (Figure  23)  the  variation  increases  up  to  about 
8  percent.   This  again  is  attributed  to  the  waves  present  in  the  flow 
field  when  the  Mach  2.0  nozzle  block  is  utilized.   The  same  trend  may 
be  observed  in  the  helium  tests  (Series  II).   With  a  Mach  1.3  primary 
stream,  static  pressure  variations  are  on  the  order  of  -   A  percent. 
With  a  Mach  2.0  primary  stream,  pressure  variations  of  up  to  9  percent 
at  about  x/a  =  5.5  (Figure  26)  may  be  seen.' 

Even  though  there  is  some  variation  of  the  static  pressure  in 
the  lateral  direction,  this  variation  is  due  to  the  presence  of  weak 
waves  in  the  supersonic  flow  field,  and  the  assumption  of  negligible 
transverse  pressure  gradient  seems  to  be  justified. 

Growth  of  the  Mixing  Region 

The  growth  of  the  shear  layers  for  each  of  the  four  test 
conditions  is  presented  in  Figure  28.   There  is  always  some  uncertainty 
in  locating  the  edges  of  the  mixing  zone.   The  range  of  uncertainty 
for  each  case  is  shown  on  the  curves  in  Figure  28. 


The  range  of  uncertainty  is  the  maximum  lateral  distance  in 
which  the  velocities  change  from  5  percent  of  their  free  stream  values 
to  the  free  stream  values. 
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The  largest  growth  is  observed  in  the  Series  II-B  tests  where 
helium  with  a  velocity  of  about  3500  ft/sec  is  mixing  with  air  flowing 
at  1300  ft/sec.   The  least  growth  rate  is  observed  in  the  case  of 
argon  (U  -  1150  ft/sec)  mixing  with  air  (U  =  1350  ft/sec),  i.e., 
Series  I-B  tests.   Thus,  with  the  two  other  configurations  showing 
the  same  trend  of  increased  mixing  zone  growth  with  an  increase  in  the 
velocity  difference,  it  is  concluded  that  as  the  velocity  difference 
between  the  two  streams  increases  so  does  the  growth  rate  of  the  mixing 
region.   This  is  consistent  with  the  well-established  fact  [1]  that  the 
growth  rate  is  a  maximum  when  one  stream  exhausts  into  a  quiescent 
medium. 

Although  curve  (a)  in  Figure  28  corresponds  to  the  test 
condition  with  the  largest  mass  flux  difference,  curve  (d)  does  not 
correspond  to  the  case  of  the  smallest  mass  flux  difference.   Therefore, 
the  same  reasoning  that  holds  for  velocity  differences  does  not  hold 
for  mass  flux  differences. 

The  curves  have  been  started  at  approximately  2  slot  heights 
downstream  of  the  exit  plane  due  to  the  lack  and  uncertainty  of  data 
at  previous  locations.   It  is  also  noted  that  if  the  curves  are 
extrapolated  to  determine  the  intercept,  none  of  the  curves  pass 
through  the  origin.   This  is  attributed  to  two  possible  reasons;  the 
first  is  the  fact  that  the  growth  of  the  mixing  zone  might  be  non- 
linear in  this  region.   The  second  and  more  probable  reason  might  be 
that  this  initial  thickness  is  due  to  the  accumulation  of  boundary 
layers  on  both  sides  of  the  splitter  plate.   It  is  the  opinion  of  the 
author  that  the  above  phenomenonis  due  to  a  combination  of  the  two 
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possibilities  rather  than  due  to  solely  one.   Of  course,  experimental 
error  could  also  have  a  significant  effect. 

Correlation  of  Velocity  Data 

A  step  function  velocity  profile  as  the  initial  profile  input 
to  the  mixing  program  could  not  be  expected  to  satisfactorily  predict 
the  mixing  in  the  near  field  since  the  profiles  are  expected  to  be 
non-similar  and  there  is  the  possibility  of  a  "wakelike"  profile 
stemming  from  the  boundary  layer  accumulation  on  the  splitter  plate. 
Thus,  it  was  necessary  to  input  measured  profiles  rather  than  assumed 
or  calculated  profiles;  and  this  was  accomplished  by  starting  the 
mixing  program  with  experimental  profiles  at  x/a  =0.57.   At  this 
axial  location,  no  "wakelike"  profiles  were  observed  in  any  of  the  test 
conditions.   As  described  in  the  Appendix,  the  input  consisted  of  velocity, 
static  temperature  and  mass  fraction  profiles  in  equal  intervals  of  Ai|i. 

The  constants  in  Equations  (A. 10),  (A. 13),  (4.14)  and  (A. 15) 
were  varied  until  the  theoretical  and  experimental  mixing  zone  growth 
rates  matched.   However,  difficulties  were  encountered  in  Prandtl's 
mixing  length  model  as  it  predicted  a  highly  non-linear  growth  rate  for 
the  mixing  zone  in  the  region  of  interest,  i.e.,  2  <_  x/a  <_  10.   No 
evident  reason  can  be  given  for  this  behavior.   It  is  speculated 
that  the  dynamic  eddy  viscosity  (pe)  is  a  very  weak  function  of  the 
transverse  coordinate.   References  [35]  and  [57]  tend  to  support  this 
speculation.   The  Prandtl  model  is  a  strong  function  of  the  transverse 
coordinate  since  it  involves  the  gradient  of  the  longitudinal  velocity. 
Hence,  the  behavior  may  be  related  somewhat  to  the  above  speculation. 
The  mixing  length  model  was  thus  eliminated  from  further  analysis. 


83 


The  remaining  three  momentum  transport  models  (Schetz,  Ferri 
and  Alpinieri)  were  correlated,  as  stated  earlier,  in  the  "very" 
near  field,  2  <_  x/a  <_  10.   The  reason  for  using  only  a  portion  of  the 
data  was  to  see  how  well  the  models  could  predict  available  mixing 
data  in  the  region  x/a  _>  10. 

Four  different  empirical  constants  were  determined  for  each  of 
the  three  models,  i.e.,  one  for  each  test  condition.   Since  the  eddy 
viscosity  is  semiempirical  in  nature,  it  is  too  much  to  expect  a  single 
correlation  to  be  valid  for  all  conditions  encountered.   These 
constants  were  then  examined  to  see  if  any  trend  could  be  observed. 
In  other  words,  it  was  desired  to  relate  the  coefficients  to  initial 
mixing  conditions.   The  only  flow  property  that  depicted  any  trend  in 
values  of  the  empirical  constants  was  found  to  be  the  ratio  of  the  mass 
flux  per  unit  area  of  the  secondary  stream  to  the  primary  stream. 
Figure  29  shows  the  values  of  the  coefficients  for  the  turbulent  eddy 
viscosity  models  as  a  function  of  the  mass  flux  ratios  of  the  two 
streams.   The  "asymptotic"  characteristic  of  coefficients  for  the  Ferri 
and  Schetz  models  as  the  mass  flux  ratio  approaches  unity  may  be 
related  to  the  fact  that  the  two  models  fail  when  the  mass  fluxes  of 
the  two  streams  are  equal.   No  reason  can  be  given  for  the  same  trend 
shown  by  the  Alpinieri  model.   It  is  also  recognized  that  more  test 
data  areneeded  in  the  region  of  unity  mass  flux  ratio  (both  less 
than  and  greater  than)  as  well  as  large  values  of  mass  flux  ratio  to 
validate  the  proposed  relation  between  the  empirical  coefficients 
and  the  mass  flux  ratios. 

The  velocity  profiles  predicted  by  each  model  as  well  as  the 
experimental  velocity  profiles  for  each  of  the  four  test  conditions 
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are  presented  in  Figures  30  through  49.   The  profiles  are  presented 
in  the  form  of  the  dimensionless  excess  velocity 


(6.1) 


and  the  dimensionless  transverse  coordinate  y/a  at  five  axial  stations. 
Only  sample  experimental  points  are  plotted  showing  the  trend  of  the 
data  to  avoid  a  cluster  of  points  in  obscuring  the  plots. 

With  all  three  of  the  eddy  viscosity  models  predicting  nearly 
the  same  growth  rate  for  the  mixing  zone,  the  Schetz  model  is 
observed  to  be  superior  in  predicting  the  velocity  profile  compared  to 
the  Ferri  and  Alpinieri  models  as  may  be  seen  in  Figures  32,  39,  43, 
47  and  49.   The  Ferri  model  tends  to  under  predict  the  velocity  profiles 
in  the  region  x/a  _>  8,  i.e.,  past  the  region  of  correlation 
(Figures  34,  38,  44  and  49).   The  Alpinieri  model  falls  in  between  the 
predictions  of  the  Ferri  and  Schetz  models.   The  better  correlation 
obtained  with  the  Schetz  model  may  be  explained  by  the  fact  that  this 
model  takes  into  account  the  velocity  and  density  profiles  in  the  shear 
layer.   By  integrating  these  profiles  at  each  axial  station  an 
"average"  value  for  the  eddy  viscosity  is  obtained.   On  the  other  hand, 
the  Ferri  and  Alpinieri  models  predict  a  turbulent  eddy  viscosity 
by  means  of  free  stream  properties  and  a  representative  width  for  the 
mixing  zone. 

Although  the  Alpinieri  and  Ferri  models  correlate  the  data 
presented  somewhat  satisfactorily,  certain  limitations  are  inherent 
in  their  formulations.   For  example,  in  both  cases,  the  mixing  zone 
width  b  is  based  on  velocity  difference;  thus  as  the  velocity 
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difference  between  two  streams  becomes  very  small,  the  mixing  zone 
width  becomes  undefined.   As  discussed  earlier,  the  Ferri  model  also 
fails  when  the  mass  flux  difference  between  the  two  streams  approaches 
zero.   Although  the  Schetz  model  also  fails  when  the  mass  flux  gradients 
in  the  flow  field  disappear,  the  former  problem  is  circumvented  since 
this  model  is  not  based  on  a  mixing  zone  width. 

Thus,  it  is  concluded  that  on  the  basis  of  formulation 
characteristics  and  the  satisfactory  correlation  of  the  experimental  data. 
the  Schetz  extension  of  the  Clauser  integral  model  for  the  turbulent 
momentum  transport  mechanism  is  superior  to  Prandtl's  mixing  length, 
Ferri' s  differential  mass  flux  and  Alpinieri's  momentum  flux  models 
in  terms  of  predicting  the  velocity  profiles  in  the  initial  region 
of  a  confined  half-jet.   It  should  be  added,  however,  that  the  Ferri 
and  Alpinieri  models  may  also  be  utilized  in  the  initial  region  with 
fairly  good  results.   Although  the  empirical  constant  in  all  three 
of  the  models  requires  adjusting  for  different  flow  configurations, 
it  is  anticipated  that  the  application  of  a  given  eddy  viscosity  model 
to  a  sufficient  quantity  of  data  will  produce  a  relationship  in  terms 
of  the  mass  flux  ratios  of  the  two  streams,  which  will  enable  calculation 
of  the  constant  for  initial  mixing  conditions. 

Correlation  of  Mass  Fraction  Data 

As  may  be  recalled  from  the  discussion  of  Section  4.5 
involving  the  turbulent  transport  of  mass,  the  approach  taken  in  this 
study  is  to  formulate  and  correlate  a  turbulent  momentum  transport 
mechanism  and  then  to  use  a  suitable  value  of  the  turbulent  Schmidt 
number  to  determine  the  speciesmass  fraction  profiles.   Thus,  the 
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correlation  technique  requires  the  employment  of  the  mixing  zone  growth 
to  determine  the  empirical  constant  in  the  eddy  viscosity  model  and 
then  the  employment  of  the  mass  fraction  profiles  to  obtain  the  most 
appropriate  value  of  the  turbulent  Schmidt  number. 

After  theeddy  viscosity  models  were  correlated  according  to 
the  procedure  discussed  in  the  previous  section,  turbulent  Schmidt 
numbers  of  0.5,  0.7  and  1.0  were  used  in  the  theoretical  analysis.   It 
was  observed  that  a  value  of  0.5  for  the  Schmidt  number  overpredicted 
the  mass  fraction  profiles  and  a  value  of  1.0  was  seen  to  be  too  high, 
i.e.,  underprediction.   Then  small  changes  in  the  neighborhood  of  the 
value  0.7  were  made.   It  was  seen  that  the  velocity  profiles  were 
insensitive  to  these  small  changes  in  the  value  of  the  turbulent 
Schmidt  number.   However,  the  mass  fraction  profiles  were  observed  to 
be  fairly  sensitive  to  changes  in  the  Schmidt  number  and  a  final 
value  of  0.7  was  chosen  for  the  Schmidt  number  as  yielding  the  best 
overall  results.   It  is  noted  that  similar  observations  were  made  in 
references  [5]  and  [35].   Whether  the  turbulent  Prandtl  number  or  turbu- 
lent Lewis  number  is  0.7  or  1.0  cannot  be  determined  from  the  data 
herein.   However,  to  be  consistent  with  the  assumption  of  constant  total 
temperature  throughout  the  mixing  region,  a  value  of  1.0  was  assumed 
for  the  turbulent  Prandtl  number.   This  determined  the  value  of  the 
turbulent  Lewis  number  on  the  order  of  1.4.   It  should  be  emphasized 
that  the  turbulent  Schmidt  number  is  sufficient  to  determine  the 
uniqueness  of  the  concentration  profiles  along  with  the  empirical 
constant  for  the  eddy  viscosity  model. 

The  helium  and  argon  mass  fraction  profiles  are  presented  in 
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Figures  50  through  69  for  all  of  the  configurations  tested.   It  was 
assumed  and  later  validated  in  Chapter  V  that  air  behaved  as  a  single 
species.  Thus,  knowing  the  secondary  stream  (helium  or  argon)  mass  frac- 
tion profile,  the  air  species  prof ile  may  be  obtained  through  the  relation 


i    .      =   1  -  a  /   ,  •,  .   \  (6.2) 

air        argon  (or  helium) 


Then  the  individual  profiles  of  oxygen  and  nitrogen  may  be  deduced  by 
the  method  employed  in  Chapter  V. 

All  theoretical  mass  fraction  profiles  were  obtained  with  a 
value  of  0.7  for  the  turbulent  Schmidt  number.   The  Alpinieri  and 
Schetz  models  in  conjunction  with  the  prescribed  value  of  the  Schmidt 
number  yielded  equally  good  results  in  terms  of  predicting  the  mass 
fraction  profiles  in  Series  I-A,  I-B  and  II-A  tests  (Figures  50 
through  64).   The  Schetz  model  was  observed  to  correlate  the  mass 
fraction  data  better  than  the  Alpinieri  model  for  the  Series  1I-B 
tests;  especially  for  axial  locations  greater  than  approximately  8 
slot  heights  downstream.   Here  the  Alpinieri  model  tended  to  over- 
predict  the  profiles  (Figures  68  and  69).   This  is  not  too  surprising 
since  the  same  trend  may  be  observed  in  the  excess  velocity  profiles 
(Figures  48  and  49).   The  Ferri  model  was  seen  to  underpredict  the 
mass  fraction  profiles  consistently.   This  is  similar  to  the  results 
obtained  for  the  velocity  profiles  with  the  Ferri  model.   The  under- 
prediction  of  the  species  prof iles  could  be  overcome  somewhat  by 
decreasing  the  value  of  the  Schmidt  number.   However,  the  same  value  of 
0.7  was  used  for  all  models  so  that  a  consistent  comparison  could  be 
made. 


Remarks  on  the  Similarity  of  Profiles 

Although  the  main  emphasis  of  this  study  is  in  the  initial  region 
of  mixing  where  the  profiles  are  expected  to  be  non-similar,  it  was 
of  interest  to  see  if  any  of  the  tests  yielded  similar  profiles  and  if 
so,  to  predict  the  onset  of  similarity.   For  this  purpose,  the  excess 
velocity  and  mass  fraction  profiles  were  plotted  in  terms  of  the 
conventional  similarity  variable  y/b.   No  similarity  features  were 
observed  for  the  Series  I-A,  I-B  and  II-A  tests.   However,  for  the  data 
of  Series  II-B  tests,  an  approach  to  similarity  was  detected  in  the 
mass  fraction  profiles.   The  similarity  plots  for  the  Series  II-B 
tests  are  presented  in  Figures  70  and  71.   It  may  be  seen  that  the 
velocity  plot  of  Figure  70  shows  no  semblence  of  similarity;  but  in 
Figure  71  it  may  be  seen  that  the  mass  fraction  data  for  x/a  =  8.44  and 
x/a  =  12.17  lie  almost  on  top  of  each  other.   From  this  it  is  concluded 
that  whEn  profiles  approach  similarity,  the  speciesmass  fraction 
profiles  exhibit  similarity  features  earlier  than  the  velocity  profiles. 
For  the  flow  conditions  of  Series  II-B  tests,  it  is  estimated  that 
the  mass  fraction  profiles  are  similar  at  approximately  12-15  slot 
heights  downstream. 


CHAPTER  VII 
SUMMARY  OF  RESULTS  AND  CONCLUSIONS 

An  investigation  of  the  turbulent  mixing  in  the  initial  region 
of  a  half-jet  composed  of  dissimilar  gas  streams  has  been  made.   An 
iso-energetic,  non-reacting  and  isobaric  binary  system  with  various 
velocity  and  mass  flux  ratios  was  studied  both  experimentally  and 
theoretically. 

The  flow  problem  was  formulated  as  an  initial  value  problem 
using  the  turbulent  boundary  layer  equations  in  conjunction  with 
phenomenological  models  for  the  turbulent  eddy  viscosity.   The  eddy 
diffusivity  was  obtained  form  the  eddy  viscosity  model  by  considering 
the  turbulent  Schmidt  number  as  a  parameter.   The  analytical  solution 
was  obtained  using  an  all-explicit  forward  marching  finite' difference 
scheme  with  the  stability  of  this  scheme  being  ensured  by  satisfying 
the  von  Neumann  stability  criterion. 

In  the  experimental  phase  of  the  study,  a  two-dimensional  test 
section  was  designed  and  built  to  be  used  in  conjunction  with  the  existing 
blow-down  wind  tunnel  facilities.   Quantitative  data  in  the  form  of  total 
and  static  pressure  measurements,  total  temperature  and  mass  fraction 
measurements  were  collected.  These  data  were  used  to  correlate  and  com- 
plement the  analytical  study. 

The  results  of  this  study  may  be  summarized  as  follows: 

1)   Constant  pressure  mixing  in  the  initial  region  of  a  half- 
jet  may  be  obtained  in  confined  supersonic  flows  with  careful  con- 
struction of  contoured  nozzle  blocks.   The  experimental  study  may  be 
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extended  to  fully  developed  mixing  zones  by  using  larger  transverse 
dimensions  and  longer  test  sections  than  used  in  this  investigation. 
The  experimental  results  also  apply  to  the  initial  regions  of  two- 
dimensional  symmetric  jets  for  which  data  are  found  to  be  quite  scarce 
in  literature. 

2)  Large  errors  in  concentration  measurements  taken  from 
binary  systems  may  occur;  the  probe  design,  sampling  technique  and 
local  turbulence  level  in  the  flow  field  have  a  significant  effect  upon 
the  results.   Satisfactory  results  were  obtained  with  the  probe  design 
and  sampling  technique  used  in  this  investigation. 

3)  Several  phenomenological  models  for  the  turbulent  momentum 
transport  were  investigated  on  a  consistent  basis.   It  is  concluded 
that  Schetz  modification  to  the  Clauser  integral  model  correlates  the 
velocity  data  in  the  near  field  very  well,  provided  the  empirical  con- 
stant is  properly  adjusted  for  each  flow  configuration.   There  is  no 
universal  constant  that  will  correlate  all  possible  flow  conditions 
simply  due  to  the  fact  that  eddy  viscosity  is  seiniempirical  in  nature. 

4)  In  the  analysis  of  different  models  for  eddy  viscosity, 
it  was  speculated  that  the  dynamic  turbulent  viscosity  (u)   is  a  very 
weak  function  of  the  lateral  coordinate.   This  may  be  a  reason  why  the 
mixing  length  model  does  not  correlate  the  data  presented.   This 
conclusion  supports  previous  work  done  in  fully  developed  mixing  zones. 

5)  In  correlating  the  mass  fraction  profiles,  it  was  observed 
that  a  constant  value  of  0.7  for  the  turbulent  Schmidt  number  gave  good 
results  for  all  flow  conditions  tested.   Therefore,  in  analyzing  tur- 
bulent mixing  problems  involving  streams  of  different  composition,  it  is 
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sufficient  to  model  a  turbulent  eddy  viscosity  and  then  retain  the 
turbulent  Schmidt  number  as  a  parameter  to  obtain  the  eddy  dif f usivities. 

6)  The  value  of  the  turbulent  Schmidt  number  being  less  than 
1.0  suggests  that  mass  appears  to  diffuse  more  readily  than  momentum. 
This  is  consistent  with  previous  work  reported  in  the  literature. 

7)  When  using  turbulent  momentum  transport  mechanisms  based 

on  the  mixing  zone  thickness  to  analyze  the  initial  region  of  a  half-jet, 
it  is  proposed  that  the  physical  width  of  the  mixing  zone  be  utilized 
rather  than  a  symbolic  half-radius  which  is  applicable  to  symmetric 
jets.   Satisfactory  correlation  results  were  obtained  with  the  Alpinieri 
model  with  the  proposed  definition  of  the  irixing  zone  thickness. 

8)  The  empirical  constant  in  a  given  eddy  viscosity  model  may 
be  related  to  initial  mixing  conditions  such  as  the  mass  flux  ratios  of 
the  two  streams.   If  a  relation  between  the  empirical  constant  and  the 
mass  flux  ratio  can  be  validated  by  large  quantities  of  data,  then   the 
constant  can  be  calculated  for  a  given  set  of  free  stream  conditions. 

9)  Boundary  layer  form  of  the  conservation  equations  is 
observed  to  be  valid  for  constant  pressure  confined  supersonic  mixing. 
The  consistency  of  the  boundary  layer  equations  was   also  checked  numer- 
ically.  It  was  observed  that  axial  derivatives  of  the  longitudinal 
velocity  were  2-3  orders  of  magnitude  smaller  than  the  corresponding 
transverse  derivatives. 

Large  quantities  of  data  are  needed  in  the  study  of  turbulent 
mixing  problems,  especially  in  the  initial  region.   As  stated  earlier, 
other  flow  conditions  than  the  ones  investigated  in  this  study  are 
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necessary  to  obtain  a  firm  relationship  between  the  empirical  constants 
and  the  initial  flow  variables.   This  study  has  dealt  with  iso-energetic 
streams;  studies  involving  streams  of  different  thermal  levels  are 
needed  so  that  solutions  and  transport  models  may  be  applicable  to  a 
broader  range  of  flow  conditions.   The  correlations  obtained  in  this  study 
are  for  the  initial  region.   It  would  be  interesting  to  see  if  these 
correlations  are  also  valid  for  "far  field"  calculations,  i.e.,  fully 
developed  mixing  zones.   Finally,  validity  of  turbulent  transport  models  in 
the  presence  of  strong  pressure  gradients  should  be  investigated  in  the 
initial  region. 
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Table  II.   Test  Conditions  and  Configurations. 


SERIES  I-A     SERIES  I-B 
(Argon-Air)     (Argon-Air) 


2.0 


1.3 


1.3 


1.3 


series  : 

[I- 

-A 

SERIES  II-B 

(Helium- 

-; 

Lr 

(Helium-Air) 

2.0 

1.3 

1.3 

1.3 

a  (in) 
(slot  ht. ) 


0.54 


0.54 


0.54 


0.54 


P 

op 

(psia) 


98+3 


41+1 


92+3 


41+1 


P 

OS 

(psia) 


37±0.5        46+0.5 


35+0.5 


46+0.5 


o 
(°R) 


535±5        535+5 


535+5 


535+5 


P 
(ft/sec) 


1710 


1300 


1710 


1300 


(ft/sec) 


1100 


1100 


3500 


3500 


*-r 


P 
(pU), 

TpUT 


0.644 


0.761 


0.846 


1.385 


2.040 


0.243 


2.69 


0.331 
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Table  III.   Comparison  of  Calculated  and  Measured 
Mass  Flow  Rates. 


Series  I-A     Series  I-B     Series  II-A    Series  II- 
/  m  /m  m  /m  m  /m  m  /m 


0.57  0.925  0.938  0.980  0.962 

2.54  0.905  0.936  0.991  0.931 

4.06  1.045  0.938  1.050  0.974 

5.56  0.944  0.946  1.020  0.990 

8.44  0.974  1.020 

12.17  0.918  1.042  1.320*  1.070 

17.70  0.895  0.950  

23.26  0.951  


Interaction  with  the  wall  and  strong  shock  patterns  were  observed 
on  schlieren  photographs. 


FIGURES 


99 


mTTTTTT 


a 


<X 


03 

c 

o 

u 
<D 
CO 


i_ 


100 


=>    8    x 


101 


E 


103 


^7 


J 


CO 

w   o 

«   X. 


Figure  6 


Close-Up  of  Test  Section  with  Two 
Mach  2.0  Nozzle  Blocks. 


Figure  7 


Close-Up  of  the  Air  Mach  2.0 
Nozzle  Block. 
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Figure  8 


Schlieren  Photograph  of  Mixing  Flow  with 
Mach  2.0  Air  and  Mach  1.3  Argon. 


Figure 


Schlieren  Photograph  of  Mixing  Flow  witl 
Mach  1.3  Air  and  Mach  1.3  Argon. 
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Figure  10  .   Probe  Drive  Mechanism  Mounted 
on  the  Test  Section 


Figure  11  .   Close-Up  of  the  Probe  Drive  Mech- 
anism with  the  Control  System 
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Figure  12  .   Close-Up  of  Static  Pressure  Plate 


Figure  13  .   Gas  Chromatograph  with  a  Molecular 
Sieve  5 -A  Column  in  Front 
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Figure  14  .   Total  Pressure  Probes. 


Figure  15  .   A  Sample  Bottle  with  Species  Sampling 
Rakes . 
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Figure  29  .   Empirical  Coefficients  for  the  Turbulent 
Viscosity  Models  as  a  Function  of  the 
Mass  Flux  Ratios  of  the  Two  Streams. 
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Figure  30 


'p-  "s 

Dimensionless  Excess  Velocity  Profile  at 
x/a=2.54  for  Series  I-A  (Argon-Air)  Tests. 


12  4 


Data 

Schetz  Model 
Alpinieri  Model 
Ferri  Model 


0.5 
Velocity    Ratio  , 


1.0 


0 


Figure  31  .   Diiriensionless  Excess  Velocity  Profile  at 
x/a=4.06  for  Series  I-A  (Argon-Air)  Tests, 
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Figure  32 


Dimensionless  Excess  Velocity  Profile  at 
x/a=5.56  for  Series  I-A  (Argon-Air)  Tests, 
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Figure  33  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=12.17  for  Series  I-A  (Argon-Air)  Tests, 
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Figure  34 


Dimensionless  Excess  Velocity  Profile  at 
x/a=17.7  for  Series  I-A  (Argon-Air)  Tests. 
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Figure  35  .   Dimensionless  Excess  Velocity  Profile  at 
>:/a=4.06  for  Series  I-B  (Argon-Air)  Tests, 
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Figure  36  .   Dimensionless  Excess  Velocity  Profile  at 
x/a=5.56  for  Series  I-B  (Argon-Air)  Tests. 
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Figure  37 


Dimensionless  Excess  Velocity  Profile  at 
x/a=12.17  for  Series  I-B  (Argon-Air)  Tests. 
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Figure  38  .   Dimensionless  Excess  Velocity  Profile  at 
x/a=17.7  for  Series  I-B  (Argon-Air)  Tests. 


132 


2.C 


1.  5 


0.5 


Data 

Schetz  Model 
Alpinieri  Model 
Ferri  Model 


0.5 
Velocity     Ratio  ,        (/) 


1.0 


Us 


'p- 


Figure  39 


Dinensionless  Excess  Velocity  Profile  at 
x/a=23.26  for  Series  I-B  (Argon-Air)  Tests. 
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Figure   40    .      Dimens:!  onless   Excess  Velocity  Profile  at 

x/a=2.54for   Series   II-A    (Helium-Air)    Tests. 
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Figure  41  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=4.06  for  Series  II-A  (Helium-Air)  Tests, 
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Figure  42  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=5.56  for  Series  II-A  (Helium-Air)  Tests. 


136 


2.0 


O   Data 

Schetz  Model 


Alpinieri  Model 

Ferri  Model 


0.5 


1.0 


Velocity     Ratio,        (f)  :       u    -   us 


Up-    US 


Figure  43  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=8.44  for  Series  II-A  (Helium-Air)  Tests, 


137 


1.  5 


-   1.0 


0.5 


Alpinieri  Model 
Ferri  Model 


Figure  44  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=12.17  for  Series  II-A  (Helium-Air)  Tests. 
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Figure  45  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=2.54for  Series  II-B  (Helium-Air)  Tests, 
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Figure  46 


Dimensionless  Excess  Velocity  Profile  at 
x/a=4.06  for  Series  II-B  (Helium-Air)  Tests. 
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Figure  47  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=5.56  for  Series  1I-B  (Helium-Air)  Tests 


2.0 


141 


1.  5 


©   Data 

Schetz  Model 


Alpinieri  Model 

Ferri  Model 


0.5 


1.0 


Velocity  Ratio,   (h  z        U  -  Us 
Up-  Us 

Figure  £8  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=8.44  for  Series  II-B  (Helium-Air)  Tests, 
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Figure  49  .   Dimensionless  Excess  Velocity  Profile  at 

x/a=12.17  for  Series  II-B  (Helium-Air)  Tests. 
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Figure  50  .   Argon  Mass  Fraction  Profile  at  x/a=2.54 
for  Series  I-A  Tests. 
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Figure  51  .   Argon  Mass  Fraction  Profile  at  x/a=4.06  for 
Series  I-A  Tests. 
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Figure  52  .   Argon  Mass  Fraction  Profile  at  x/a=5.56  for 
Series  I-A  Tests. 
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Figure  53  .   Argon  Mass  Fraction  Profile  at  x/a=12.17  for 
Series  I-A  Tests. 
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Figure  54  .   Argon  Mass  Fraction  Profile  at  x/a=17.7  for 
Series  I -A  Tests. 
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Figure  55  .   Argon  Mass  Fraction  Profile  at  x/a=A.06for 
Series  1-B  Tests. 
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Figure  56  .   Argon  Mass  Fraction  Profile  at  x/a=5.56  for 
Series  I-B  Tests. 
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Figure  57  .   Argon  Mass  Fraction  Profile  at  x/a=12.17  for 
Series  I-B  Tests. 
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Figure  58  .   Argon  Mass  Fraction  Profile  at  x/a=17.7  for 
Series  I-B  Tests. 
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Figure  59  .   Argon  Mass  Fraction  Profile  at  x/a=23.26  for 
Series  I-B  Tests. 
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Figure   60    .      Helium  Mass  Fraction  Profile   at  x/a=2.54for 
Series   II-A  Tests. 
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Figure  61  .   Helium  Mass  Fraction  Profile  at  x/a=4.06  fcr 
Series  II-A  Tests. 


155 


1  .5 


0.5 


A   Data 

Schetz  Model 

Alpinieri  Model 
Ferri  Model 


I J L. 

0.5 
Mass     Fraction  ,  Q( 


1.0 


Helium 


Figure  62  .   Helium  Mass  Fraction  Profile  at  x/a=5.56  for 
Series  II-A  Tests. 
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Figure  63  .   Helium  Mass  Fraction  Profile  at  x/a=8.44  foi 
Series  II-A  Tests. 


157 


2.0 


1 . 0 


A   Data 

Schetz  Model 
Alpinieri  Model 
Ferri  Model 


0.5 
Mass  Fraction  ,  Q( 


Figure  64  .   Helium  Mass  Fraction  Profile  at  x/a=12.17  for 
Series  IT-A  Tests. 
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Figure  65  .   Helium  Mass  Fraction  Profile  at  x/a=2.54  for 
Series  II-B  Tests. 
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Figure  66  .   Helium  Mass  Fraction  Profile  at  x/ac4.06  for 
Series  II-B  Tests. 
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Figure  67  .   Helium  Mass  Fraction  Profile  at  x/a=5.56  for 
Series  II-B  Tests. 


161 


[ 


1.5 


0.5  - 


A   Data 

Sche.tz  Model 

Alpinieri  Model 

Ferri  Model 


J L 


J L 


0.5 

Mass     Fraction  ,  (Y 

Helit 


-J 
1.0 


Figure  68  .   Helium  Mass  Fraction  Profile  at  x/a-8 .44  f or 
Series  II- B  Tests. 
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Figure  6y  .   Helium  Mass  Fraction  Profile  at  x/a=12.17  for 
Scries  II-B  Tests. 
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Figure  70  .   Velocity  Similarity  Plot  for  Series  II-B  Tests 
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Fipurr.  71.  Mass  Fraction  Similarity  Plot  for  Series  II-B  Tests 
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APPENDIX 

Description  of  the  Program 

The  computer  program  consists  of  five  routines  for  which 
the  flow  diagrams  are  shown  in  Figures  72  through  76.   The  function 
of  each  routine  may  be  summarized  as  follows: 

1)  The  MAIN  routine  is  used  for  input,  defining  constants, 
initialization  2nd  to  link  the  rest  of  the  subroutines  together. 

2)  Subroutine  PRELIM  computes  the  parameters  necessary 
to  implant  the  finite  difference  scheme  and  also  controls  output. 

3)  Control  passes  from  subroutine  PRELIM  to  subroutine  SHEAR 
which  computes  the  thickness  of  the  mixing  zone.   Then  multiple 
entries  are  utilized  to  enter  the  same  subroutine  depending  on  the 
model  for  the  turbulent  eddy  viscosity. 

4)  Once  the  parameters  and  viscosity  are  computed,  control 
is  passed  to  subroutine  SOLVE  which  solves  the  explicit  finite 
difference  equations  (3.41)  through  (3.43)  and  (3.45)  through  (3.47). 

5)  After  a  step  in  the  streamwise  direction  is  taken  by 
solving  the  difference  equations,  subroutine  SUPP  is  entered.   SUPP 
is  a  complementary  subroutine  to  SOLVE.   It  takes  the  values  of 
calculated  total  enthalpy,  velocity  and  mass  fractions  from  SOLVE  and 
computes  static  mixture  and  species  enthalpies ,  mixture  specific 
heats,  molecular  weights,  densities  and  temperatures  at  each  new 
grid  point. 

The  program  terminates  and  calls  for  new  initial  profiles 
when  it  has  "marched"  to  the  specified  streamwise  location,  XMAX. 
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Input  to  the  program  is  a  finite  number  (<_   100)  of  grid 
points  in  the  lateral  to  flow  (v0  direction.   The  input  grid  points 
are  initial  velocity,  static  temperature  and  specie  mass  fraction 
profiles  in  equal  intervals  of  AiJ;. 

All  major  variables  used  in  the  program  are  defined  at  the 
beginning  of  the  Fortran  listing.   It  is  hoped  that  these  definitions 
in  conjunction  with  numerous  comment  cards  in  the  program  enable  the 
reader  to  follow  the  Fortran  program  with  relative  ease. 

As  the  program  executes  the  computations,  additional  points 
in  the  lateral  direction  are  added  one  at  a  time  to  the  grid  whenever 
it  is  necessary.   The  criterion  for  adding  points  to  the  mesh  is 


"n.MPSI  "  "n.MPSI-1 
^n.MPSI 


>  0.001  (A.l) 


where  0,   denotes  enthalpy,  velocity  or  the  element  mass  fractions  and 
MPS I  is  the  total  number  of  ty   grid  points. 

The  points  are  added  until  (2k-l)  points  (where  k  is  the  number 
of  points  originally  input)  are  obtained  in  the  grid;  at  which  time 
the  4/   grid  increment  is  doubled  and  the  number  of  points  are  reduced 

to  k,  i.e.,  only  \p   point  numbers  1,3,5,7 ,k  are  retained  by  the 

program.   All  ty   grid  points  not  in  use  take  on  primary  stream  values. 

Physical  Coordinates 

Once  a  step  in  the  streamwise  direction  is  taken,  the  physical 
coordinate  corresponding  to  all  4>   points  is   calculated  through  the 
relations 
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+  Ml 

'm+1   ym   2 


1   +    1 


(pu)     (pu) 

m       m-1 


(A. 2) 


Step  Size 

As  can  be  seen  from  Equation  (3.53)  there  is  no  restriction 
on  the  grid  size  AiJ/.   However,  A'^  has  to  be  as  small  as  feasible  to 
ensure  finite  difference  consistency  and  accuracy,  i.e.,  C  in  Equation 
(3.^0)  has  to  be  small.   An  approximate  relation  between  Aiji  and  the 
number  of  grid  points  may  be  explained  by  the  following  analysis: 

The  von  Mises  transformation  is 

3d 

r1-  =  pu 
oy 

y       ^ 

Je  re 

pudy  =     dip 
0         0 


Hence, 


-\    » 


>udy 
0 

where  the  subscript  e  refers  to  the  lateral  coordinate  at  which  flow 

consists  of  only  the  primary  stream. 

Now  as  an  approximation  a  step  function  velocity  profile  as 

the  initial  profile  at  x  =  0  is  assumed.   Then,  pure  secondary  gas  (argon 

or  helium)  exists  between  y  =  0  and  y  =  a  (where  a  is  the  exit  height 

of  the  secondary  stream,  i.e.,  slot  height).   Hence,  pu  is  constant 

from  y  =  0  to  y  =  a.   Therefore 

ij>  =  (pu)  a 
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The  external  streamline  ty     may  be  related  to  the  number  of 
input  grid  points  by 

(k  -  l)A<{i  =  i{»e  (A.  3) 

a(pu) 
Then  A^  -  k  _  ±  (A. A) 

The  number  of  input  grid  points  can  thus  be  estimated  for 
a  given  value  of  A^.   Obviously,  for  a  given  configuration,  the 
smaller  the  step  size  A^,  the  larger  will  be  the  number  of  grid 
points.   This  would  result  in  an  inefficient,  slow  and  also  expensive 
program. 

For  all  cases  under  study,  a  lateral  step  size  of  5x10   was 
used.   This  resulted  in  the  use  of  approximately  20  grid  points  for 
the  helium  tests  and  50  grid  points  for  the  argon  tests  with  about 
5  initial  points  in  the  mixing  zone. 

It  should  be  noted  that  Equation  (A. 4)  is  just  an  approximate 
analysis  to  determine  an  appropriate  value  A^.   After  a  value  for 
L\p   is  chosen,  the  actual  number  of  grid  points  is  determined  by  the 
experimental  curves  of  p,  U,  T  and  a..      Inside  the  mixing  region  equal 
increments  of  \p  do  not  yield  equal  increments  of  the  physical  coordinate 
y.   Thus,  some  iteration  with  Equation  (A. 2)  is  required  when  putting 

in  the  actual  initial  profiles  at  some,  x  =  x  . 

-4 
The  streamwise  step  size  used  was  Ax  =  4.5  x  10   feet. 

However,  the  stability  requirement  of  Equation  (3.53)  imposes  a 

restriction  on  Ax.   Thus,  given  a  value  of  Ai|>  (i.e.,  0.005),  the 

upper  limit  on  Ax  was  calculated  for  each  grid  point.   The  minimum 

-4 
of  the  calculate  Ax  was  compared  with  the  value  of  4.5  x  10   feet. 
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The  smaller  of  the  two  values  was  chosen  for  the  streamwise  step 
size. 

One  set  of  helium  and  argon  data  were  run  with  streamwise 

-3  -5 

step  sizes  of  5  x  10   and  5  x  10   feet  to  observe  the  effect  of 

C,  i.e.,  Equation  (3.40),  on  the  results.   No  differences  were 

-4 
detected  in  either  case  in  comparison  to  a  step  size  of  4.5  x  10 

feet.   It  is  deemed  that  with  the  use  of  the  stated  step  size, 

is  insignificant  in  the  streamwise  region  of  interest,  i.e.,  x/a  <_  20. 

Enthalpy  Fits 

To  initially  implement  the  finite  difference  scheme,  total 
enthalpies  had  to  be  calculated  from  the  input  data  of  temperature, 
velocity  and  mass  fractions.   Then,  once  the  computations  had  started, 
temperatures  had  to  be  determined  from  the  newly  calculated  values  of 
total  enthalpies. 

The  species  enthalpies  (Btu/lbm)  were  curve  fitted  by  means  of 
the  following  equation: 

h.  =  A,  +  B.T       T.    <  T  <  T,  .  (A. 5) 

ill        low  —   —  nigh 

T  being  the  absolute  temperature  in  °R,  T  ~   and  T,  .  ,  being  the 
range  of  validity  of  the  curve  fit.   The  value  of  the  constants  for 
each  species  are  given  in  the  following  table. 

These  fits  are  of  data  from  reference  [86]  and  are  based  on 
a  reference  temperature  of  400°R. 

It  is  also  noted  that  the  constants  B.  correspond  to  the 


specific  heat  at  constant  pressure  of  species i. 
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Table  IV.   Enthalpy  Fit  Coefficients 


Max.  Error 


Soecies 

A. 
l 

B. 
l 

Tl 

(°R) 
ow 

250 

Th: 

Lgh 
550 

Curve  Fit 

i   o2 

-0.47 

0.2175 

0.05% 

2    N2 

-0.11 

0.2481 

250 

550 

0.05% 

3   A 

-0.05 

0.1244 

325 

550 

0.1% 

4   He 

-0.50 

1.2413 

325 

550 

0.1% 
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R.ead  Input 


Compute  Mixture  Molecular 
Weights,  Densities  and  Specific  Heats 


Compute  Initial  Enthalpies 
(Total,  Static,  Species) 


Compute  ijj  and  Fill  Unused 
Grid  Points  with  Free  Stream 
Values 


Output  Initial  Values 
I 


Check  Input 


Terminate  if  Negative 
P,  T  or  a. 


Call  Subroutines  and  Loop 


Figure  72  .   Flow  Diagram  for  MAIN  Routine, 
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Compute  Physical 
Coordinates 


Call  SHEAR 


Compute  puE 


Determine 

Whether  to 

Output 


NO 


Return- 


NO 


-  Normal  Output 

r 

Determine 

Whether  Intermediate 

Dump  is  Requested 

YES 

Intermediate  Dump 


Punched  c  YES 
Output 


Determine 

Whether  Punched  * 

Output  is  Requested 


Figure  73  .   Flow  Diagram  for  Routine  PRELIM. 
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Compute  Mixing 
Zone  Thickness 


Branch  to  Appropriate 
Eddy  Viscosity  Model 


Prandtl's  Mixing 
Length  Model 


Schetz's  Integral 

Model 


Ferri's  Differential 
lss  Flux  Model 


Alpinierie  s 
•turn  £  Mass 
Flux  Model 


v 
Return 


Figure  74  .   Flow  Diagram  for  Routine  SHEAR. 
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Determine  Step  Size 


Solve  Explicit  Equations 
at  the  Wall 


Determine  Whether  Both 

Prandtl  &  Lewis  Numbers 

are  Unity 


Solve  Interior  Grid 
Points 


Call  SUPP 


Determine  Whether  Negative 
Temperature  or  Species  is 
Obtained 


NO 


Determine  Whether 


YES 


Return  With 
Termination 
Command 


->  Halve  Mesh 


to  naive 

nesn 

YES 

NO 

\ 

Test  for  Zero  Slope 

at  Edge 

YES 

NO 

Add  Mesh  Points 
1 

1 
Return  -g   

s 

Figure  75  .   Flow  Diagram  for  Routine  SOLVE. 
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Compute  Mixture  Static 
Enthalpy 


Compute  Temperature 


Compute  Mixture  Molecular 
Weights  and  Specific  Heats 


■Compute  Density  and 
Specie  Enthalpies 


Figure  76  .   Flow  Diagram  for  Routine  SUPP. 
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DEFINITION    OF   VARIABLES   IN    COMPUTER  PROGRAM 


{DATA   SET    ID   CARD-ANY   72   CHARACTERS 
(DETERMINES   TUR3ULENT   TRANSPORT  MECHANISM  MODEL 

(=1   PRANDTL   I) 

(=2  MODIFIED  FERRI) 

(=4   SCHET2) 

(=5   MODIFIED  ALPLNIERI) 
INUM3ER   OF   GRID   POINTS    IN    PSI   DIRECTION 
j CONTROLS   PRINTED   OUTPUT 

(=0  NORMAL   OUTPUT) 

(/O   NORMAL   OUTPUT   +    INTERMEDIATE  DUMP) 
INUM3ER   OF   SPECIES 
iCONTROLS    PUNCHED    OUTPUT 

(=0   NO   PUNCHED    OUTPUT) 

(/O    PUNCHED   OUTPUT) 
: CONTROLS  WHETHER   TO   OUTPUT   IN    INTERVALS 

(=0    OUTPUT   LN    INTERVALS) 

(j*0    OUTPUT  @  PRESCRIBED   X-LOCATIQNS) 
: PRESCRIBED   X-STATI0N3   FOR   OUTPUT    (FT) 
lENTHALPY   FIT  COEFFICIENT    (BTU/LBM) 
:ENTHALPY   FIT  COEFFICIENT    (BTU/lBM-R) 
j ENTHALPY   FIT  COEFFICIENT    (BTU/LBM-R**2) 
i INCREMENT    IN   FEET   PROGRAM   IS    TO   OUTPUT 
i INITIAL  VALUE   OF  X    IN   FEET 
lL  VALUE   OF  X    LN    FEET 

{prandtl  number 
: lew is  number 

i mesh  increment  in  psi  direction 
:lity  step  size  regulator 

iSTATIC    PRESSURE   LN    L3FZfT**2 
j SLOT  HEIGHT    (FT) 

lEMPIRICAL  TURBULENT  VISCOSITY   COEFFICIENTS 
i STATIC    TEMPERATURE    (R) 
{VELOCITY    (FT/SEC) 

iN ON-DIMENSIONAL  VELOCITY   RATIO    (U-UJ)/(UE-UJ) 
j  SPEC  IS   MASS   FRACTIONS 
{MOLECULAR  WEIGHTS    OF   SPECIES 
{MOLECULAR  WEIGHT   OF   MIXTURE 
{MIXTURE   DENSITY    (SLUG3/FT**3 ) 
{MIXTURE   TOTAL   ENTHALPY    (BTU/LBM) 
{MIXTURE   STATIC    ENTHALPY    (BTU/LBM) 
) {SPECIES    STATIC    ENTHALPY    (BTU/LBM ) 
i(U*U)/2 
{STREAM  FUNCTION    (SLUGS/FT-SEC ) 

ffiER    OF   STEPS    m    X-DIRECTI 
{DELTA  X    INCREMENT    (FT) 


c 

DEFINITION    ( 

c 

c 

TITLE (18) 

c 

MODEL 

c 

c 

c 

c 

c 

MPS  I 

Q 

IN  DUMP 

c 

c 

c 

NOSPSC 

c 

N PUNCH 

c 

c 

IOUT 

c 

c 

c 

XBAR(9) 

c 

kW 

c 

BW 

Q 

C(4) 

c 

Y  p  tj  7^  j  rn 

c 

XINIT 

r 

XFINAL 

c 

PRNDTL 

c 

XLEWIS 

c 

DSLPSI 

c 

XMPS 

c 

P 

c 

HEIGHT 

Q 

CI.,.. .06 

c 

T(200) 

c 

U(200) 

c 

PHI (200) 

c 

ALFA (^,2 00) 

c 

WTMOL(^) 

c 

WTMIX(200) 

c 

RH0(200) 

c 

HT0T(230) 

c 

HSTAT(200) 

c 

HSPEC(k,200 
U3Q(200) 

c 

c 

P3I(200) 

c 

ISTEP 

c 

DX 
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C  JEXIT  i CONTROLS   TERMINATION    OF   PROGRAM 

C  Y(200)  iY   DISTANCE   FROM  LOWER  WALL    (FT) 

C  E?SI(200)         tTURBULENT  VISCOSITY    (L3F-3EC/FT**2) 

C  SG(200)  iRHO*U*EPSI 

C  CP3AR(200)       : MIXTURE   SPECIFIC   HEATS    (BTU/L3M) 

C  SIGMA (200)       : CHECK  SUM  FOR  MASS   FRACTIONS 

C  WU(200)  : VELOCITY    (WORKING  ARRAY) 

C  WALFA(^,200) J MASS   FRACTIONS    (WORKING  ARRAY) 

C  WH(200)  j TOTAL  ENTHALPY    (WORKING  ARRAY) 

C 

C 

DIMENSION  TITLE(l8),A(M,3(^),C(iOtT(200),U(200). 
1WTM0L(^),WTMIX(200),RH0(200),HT0T(200),H3TAT(200), 
2USQ(200),PSI(200),Y(200),E?3I(200),SG(200),C?3AR(200), 
3WU(200),WALFA(4,200),WH(200),ALFA(^,200), 
^HSPECd,  200), SIGMA  (200) 
C 

C  OMMCN  TITLE ,  MPSI ,  INDUMP  ,N  OSPEC  ,NPUNCH ,  A , 3 ,  C  ,  XPRIN T ,  XIN  IT 

COMMON  XFINAL,PRNDTL,XLEWIS,DELPSI,XMPS,P,C1,C2,C3,C^,T,U 

COMMON  ALFA , WTMOL ,WTMIX ,  RHO , HTOT ,HSTAT ,HSPEC , USQ , PSI , ISTSP 

COMMON  PCNT,DX,JEXITfY,EPSI,SG,CPBAR,SIGMA,WU,WAI»FA,WH,NPSI 

C  OMMON    X ,  MIN I T ,  MHALF ,  IPAGE 

COMMON    X3AR(9 ) , IX3AR,PHI(200 ) ,BZ,C5,C6, MODEL, HEIGHT 
C 
C 

C  INPUT   DATA 

C 
100        READ    (5,i,END=99)       (TITLE( I) ,1=1,15) 

READ (5 i 1969)    MODEL 

READ (5 , 2 )    MPSI , INDUMP ,N OSPEC ,NPUNCH, I OUT 

I?    (MPSI.EQ.O)    GO   TO  99 

IF    (IOUT.NE.O)    READ    (5,1970)      X3AR 

DO   k   1=1 ,N OSPEC 

READ    (5,3)   A(I),3(I),C(I) 
k  CONTINUE 

READ    (5,5)       ( WTMOL ( I ),I=1,N OSPEC) 

READ    (5,6)   XPRINT,XINIT,XFINAL 

READ    (5,7)    PRNDTL,XLEWIS,DELPSIfX:.PS 

READ    (5,8)   P, HEIGHT 

READ    (5,81)    C1,C2,C3,C^,C5,C6 

READ    (5,9)    (T( I), 1=1, MPSI) 

READ    (5,9)    (U( I), 1=1, MPSI) 

DO   12    1=1, MPSI 

READ    (5,11)    (ALFA(J,I),J=1,N0S?SC) 

12        continue 

c 

C  COMPUTE  MIXTURE   MOLECULAR  WEIGHTS  AND  DENSITIES 

C 

DO   14   1=1, MPSI 

SUM=0.'0 
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DO    13    J=l,NOSPEC 

SUM=S  UM+ALFA  ( J ,  I  )/WTMOL  ( J  ) 
13  CONTINUE 

WTMIX(I)=1./SUM 

RH0(I)=(P*WTMIX(I)A9719.1*T(I)) 
\k         CONTINUE 

!  COMPUTE  ENTHALPIES 

DO   16    I=1,MPSI 

SUM=0.0 

DO   15   J=l,NOSPEC 

HSPEC(J,I)=A(J)+(3(J)+C(J)*T(I))*T(I) 

SUM=SUM+H3PEC  ( J ,  I )  *ALFA  ( J ,  I ) 

15  CONTINUE 
HSTAT(I)=3UM 
USQ(I)=U(I)*U(l)/50073. 
HTOT ( I )=HSTAT ( I )+USQ ( I ) 

16  CONTLNUE 

COMPUTE  MIXTURE  SPECIFIC  HEATS 

DO  30   1=1, MPS I 
PHI(I)=(U(I)-U(1))/(U(MPSI)-U(1)) 

SUM=0.0 

DO   31    J=1,N0SPSC 

IF    (C(J))    ^0,M,^0 
M  SUBfeSUM+B  ( J )  *ALFA  ( J ,  I ) 

GO   TO  31 
kO  SUM=SUM+(B(J)+2.*C(J)*T(l))*ALFA(J,I) 

CONTINUE 

CPBAP-SUM 
30  CONTINUE 

COMPUTE  PSI   AND  FILL  UNUSED   GRID  POINTS 
3  WITH  FREE  STREAM  VALUES 

DO   17    1=1,200 

XI=I-1 

PSI=XI*DELPSI 

17  CONTINUE 

DO   18    1= MPS I, 200 

RH0(I)=RH0(MPSI) 

T(I)=T(MPSI) 

HTOT(I)=HTOT(MPSI) 

U(I)=U(MPSI) 

PHI(I)=PHI(MPSI) 

CPBAR ( I ) =C P3AR ( MPS I ) 

WTMIX ( I )=WTMIX(MPSI ) 

USQ(I)=USQ(MPSI) 

DO   18    J=1,N0SPEC 
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ALFA  (J,  I  )=ALFA  (J ,  MPSI ) 
13  CONTINUE 

ISTEP=0 

NPSI=MPSI-1 

JEXIT=0 

X=XINIT 

DX=0.0 

I?AGS=0 

IX3AR=1 
C 

C  OUTPUT   0?   INPUT 

C 

WRITS    (6,19) 

WRITS    (6,1)    (TITLE(I), 1=1,15) 

WRITE    (6,1969)    MODEL 

WRITS    (6,2)    MPSI,INDUMP,NOS?SC,NPUNCH,IQUT 

IF    (IOUT.NE.O)   WRITE   (6,1970)    X3AR 

DO  20    I=1,N0SPSC 

WRITS    (6,3)   A(I),3(I),C(I) 
20  CONTINUE 

WRITS    (6,5)    (WTM0L(I),I=1,N0S?SC) 

WRITE    (6,6)    X?RINT,XINIT,XFINAL 

WRITS    (6,7)    PRNDTL,XLF,VI3,D3LP3I,XM?S 

WRITS    (6,8)    P, HEIGHT 

WRITS    (6,31)    C1,G2,G3,G^,C5,G6 

WRITS    (6,9)    (T(l), 1=1, MPSI) 

WRITE   (6,9)    (U(I),I=1,MPSI) 

DO   22    1=1, MPSI 

WRITS    (6,11)     (ALFA(J,I),J=1,N0S?SG) 
22  CONTINUE 

C 

C  CHECK   INPUT 

C 
29  IF    (P.LE.O.O)      GO    TO   2h 

DO   23    1=1,  LIPS  I 

IF    (T(I).LE.O.O)       GO   TO   2k 

DO   23    J=1,N0SPEC 

IF    (ALFA(J,I).LT.0.0.0R.ALFA(J,I).GT.1.0)      GO   TO  2k 


23 

CCNTITiUS 
GO  TO  26 

24 

WRITS  (6,25) 
GO  TO  100 

26 

CALL  PRELIM 

IF  (JSXIT.GT.l)   GO 

TO  100 

CALL  SOLVE 

GO  TO  26 

C 

; 

c 

1 

FOR; AT  (13A4-) 

2 

FORMAT  (18,13,^110) 
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3  FORMAT  (3F10.5) 

5  EORMAT  (J4.F1 0.3) 

6  FORMAT  (3F10.5) 

7  FORMAT  (2F10.,2,F10.5.F10.1) 

8  FORMAT  ( Fl 0 . 1 , 6F1 0 . 5 ) 
81  FORMAT  (6F10.40 

9  FORMAT  (T10,7F9.1) 

1969  FORIvIAT  (T10,I1) 

1970  FORMAT  (9F8.5) 

11  FORMAT  (T10.4F10.4) 

19  FORMAT  (•!•  ) 

25  FORMAT  (•O1 ,//T21, ******  PROGRAM  TERMINATED  DUE  TO  ERROR' 

1,'    IN    P,T,ALFA   ->****»  ) 

99  STOP 
END 


SUBROUTINE  PRELIM 

THIS   SUBROUTINE  COMPUTES   THE  PARAMETERS  NECESSARY   TO 
IMPLANT   THE  FINITE  DIFFERENCE  SCHEME.      THE   OUTPUT   IS  ALSO 
CONTROLLED  BY   PRELIM 

DIMENSION   TITLS(18),A(^),3(^)SC(^),T(200)SU(200), 
1V/TMOL  ( k )  ,WTMIX  ( 200  ) ,  RHO  ( 200  )  ,HTOT  ( 200  ) ,  HS TAT  ( 200  ) , 
2USQ(200),PSI(200),Y(200),EPSI(200),SG(200),CPBAR(200)I 
3WU ( 200 ) ,WALFA (4, 200 ) ,WH (200 ) , ALFA (4, 200 ) , 
iffiSPEC(4t200),SIGMA(200) 

C  OMMON    TITLE ,  MPS I ,  IN  DUMP  ,N  CSPEC  ,N  PUN  CH ,  A ,  B ,  C  ,  XPRIMT ,  X  IN  IT 

COMMON   XFINAL,PRNDTL,XLEWIS,DELPSIlXMPS,P>ClfC2,C3,C4fTfU 

C OMMON   ALFA , WTMOL , WTMIX , RHO ,HTOT , HSTAT, HSPEC , USQ , PS I , ISTSP 

C OMMON    PCNT, DX , JSXIT , Y , SPSI , SC- , CP3AR ,SIGMA ,WU,WALFA ,WH ,NPS I 

C  OMMON   X , MIN  IT , MHALF , IPAGE 

COMMON    ?HI(200) ,3Z,C5,C6, MODEL, HEIGHT 

C  OMMON    XBAR ( 9 ) , IX3AR 

COMPUTE  PHYSICAL   COORDINATES 

Y(1)=0.0 

DO   1000    1=2, MPS I 

Y(I)=Y(I-l)+DELPSI*(l./RH0(I)/U(l)+l./RH0(l-l)/U(I-l))/2. 

CONTINUE 

COMPUTE  TURBULENT  VISCOSITY 

CALL  SHEAR 

IF    (MODEL. SQ.l)      CALL  SHEAR1    (&19W 


c 
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IF  (MODEL. EQ. 3)  CALL  SHEAR3  (&1984) 
IF  ( MODEL. EQ.*0  GALL  SHEARS  (&1934) 
IF    (MODEL. EQ. 5)    CALL  SHEAR5    (&1984) 


198^      ISTEP=I3TEP+1 

DO  2000  1=1, MPS I 

SG ( I )=RHO ( I ) *U ( I ) *E?S 1(1) 
2000  CONTINUE 

Q 

C  DETERMINE  WHETHER  TO   OUTPUT 

G 

IF    (JEXIT)      300,100,300 
300        CONTINUE 

IF    (X.GS.XFINAL)    GO   TO   101 

IF    (IOUT.NE.O)    GO   TO   102 

IF    (XPRINT.LE.PCNT)    GO  TO   100 

PCNT=PCNT+DX 

GO  TO   350 
102        IF    (X.LT.XBAR(IXBAR))    GO  TO   350 

IXBAR=IXBAR+1 

GO  TO   100 
101        JEXIT=2 

NN=1 
60         WRITE   (6,55) 

KK=1 

WRITE    (6,62) 

IF    (N0SPEC.EQ.3)    GO   TO  65 

WRITS    (6,63) 
65  CONTINUE 

DO  70    I=NN,MPSI 

WRITE    {6,67)    I,HSTAT(I),(HSPEC(J,I),J=1,N0SPEC) 

KK=KK+1 

IF    (KK.LT.25)    GO   TO  70 

NN=I+1 

GO   TO   60 
70  CONTINUE 

C 

C  COMPUTE  SIGMA  ALFA  AS  A   CHECK  SUM 

C 
100        DO   80    I=1,MPSI 

SIGMA(I)=0.0 

DO  80   J=1,N0SPEC 

S IGMA ( I ) =S IGMA ( I ) +ALFA ( J , I ) 

60  continue 
pc:;t=o.o 
ipage=ipage+1 


OUTPUT 

::n=i 

WRITE  (6,210)  (TITLS(I),I=1_,15) 
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IF  (MODEL. EQ.l)  WRITE  (6,1971) 

IF  (MODEL.EQ.3)  WRITE  (6,1973) 

IF  ( MODEL. EQ. if)  WRITS  (6,197*0 

IF  (MODEL. EQ. 5)  WRITS  (6,1975) 
WRITS  (6,211) 

xdum=x/heigkt 

dxdum=dx/height 

bzdum=bz/height 

writs  (6,212)  xdum,pf?rndtl,xlewi3,3zdum,dxdum,istep,ipag] 

WRITS  (6,5212) 
216   WRITE  (6,213) 

KK=1 

DO  21 4  I=NN,MPSI 

YDUM=Y(l)/nEIGHT 

WRITE    (6,215)    I,P3I(I),U(I),T(I),RH0(I),EPSI(I),PHI(I), 
1YDUM.I 

KK=KX+1 

IF    (KK.LT.25)    GO   TO    21^ 

NN=I+1 

WRITS    (6,5211) 

GO   TO   216 
214-        CONTINUE 

NN=1 

WRITS    (6,210)     (TITLE(I),I=1,15) 

IF    (MODEL. EQ.l)   WRITE    (6,1971) 

IF    (MODEL. EQ. 3)    WRITE    (6,1973) 

I?   ( MODEL.. EQ.40   WRITS    (6,197^) 

IF    (MODEL. EQ. 5)   WRITS    (6,1975) 

WRITE    (6,310)   XDUM.IPAGE 

WRITS    (6,5212) 
325       WRITE   (6,311) 

IF    (N0SPEC.LE.3)    GO   TO  315 

WRITE    (6,312) 
315         CONTINUE 

KK=1 

DO   32^  I=NN,MPSI 

WRITS    (6,320)    I,SIGMA(l),HTOT(I),CPBAR(l),WTMIX(l), 
1(ALFA(J,I),J=1,N0SPEC) 

KK=KK+1 

IF    (KK.LT.25)    GO   TO   32^ 

NN=I+1 

WRITS    (6,5211) 

GO   TO   325 
32^        CONTINUE 
C 

C  DETERMINE  WHETHER   INTERMEDIATE  DUMP   IS   REQUIRED 

C 


IF    (INDUMP.EQ.O)    GO   TO   330 
NN=1 
3^1        WRITE    (6,210)     (TITLE (I), 1=1, 15) 
IF    (MODEL. EQ.l)    WRITS    (6,1971) 
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IF    (MODEL. EQ. 3)    WRITE    (6,1973) 

IF    ( MODEL.  EQ.Jf)    WRITE    (6,197*0 

IF   (MODEL.EQ.5)  WRITE   (6,1975) 

WRITS    (6,62) 

IF    (NOSPEC.LE.3)    GO  TO  335 

WRITS    (6,63) 
335        CONTINUE 

KK=1 

DO   3^0    I=NN,MPSI 

WRITS    {6t6?)    I,HSTAT(l),(HS?EC(J,I),J=i,NC3PSC) 

KK=KK+1 

IF    (XX.LT.25)    GO   TO  3^0 

NN=I+1 

GO   TO   3 ij-1 
3^0        CONTINUE 
G 
C  DETERMINE  WHETHER   PUNCHED   OUTPUT    IS   REQUIRED 

"330        IF    (NPUNCH.EQ.O)    GO   TO   350 

12=1 

KK=1 

L=l 

JJ=7 
353        I?    (JJ.LT.IJPSI)    GO  TO   351 

JJ=MPSI 
351        WRITS    (7,352)    IZ,ISTEP,XK,(T(LL),LL=L,JJ) 

IF    (JJ.EQ.MPSI)    GO   TO   ^50 

L=JJ+1 

JJ=JJ+7 

KK=KK+1 

GO  TO  353 
MO        CONTINUE 

IZ=2 

KK=1 

L=l 

JJ=7 
^-52         IF    (JJ.LT.MPSI)    GO   TO   if5l 

JJ=MPSI 
k-51        WRITS  (7,352)  IZ,I3TEP,XX,(U(LL),LL=L,JJ) 

IF  (JJ.EQ.MPSI)  GO  TO  550 

L=JJ+1 

JJ=JJ+7 

KX=XX+1 

GO  TO  k-52 
550   CONTINUE 

IZ=3 

DO  552  1=1, MPS I 

KK=I 

WRITS  (7,551)  IZ,I3TEP,XX,(ALFA(J,I), J=1,N0SPEC) 
552   CONTINUE 
C 
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C 
G 
55    FORMAT  ('1' ,T9, ******  DUMP  OF  MIXTURE  AND  SPECIES  STATIC 

1*  ENTHALPIES  AT  LAST  STEP  DUE  TO  NORMAL  PROGRAM', 

2'  TERMINATION  ******  ) 

62  FORMAT  (//T2,'PT.» ,T10,' STATIC  ENTHALPY'  ,T36, • HALPHA (1 ) • , 
1T55 , ' HALPHA(2) ' ,T76 , 'HALPHA (3 )  ■  ) 

63  FORMAT  ( •  +  •  ,  T96 , '  HALPHA  W) 
67           FORMAT  ("O* fT2,i3tlP5220.6) 

210  FORMAT  ( • 1 • , T25  * 1 3Xk ) 

211  FORMAT  ('0'  ,T1^5»X/A'  tT32, *  PRESSURE*  ,T48,  •  PRANDTL'  ,T50, 
1» LEWIS* ,T75,,BZ/A» ,T92,'STEP  SIZE' ,T107 ,* STEPS' ,T115, 
2' PAGE' ) 

212  FORMAT  ( ' 0'  1P2E20.6 ,T50,0PF^.2,T6l ,F±. 2 ,T71 ,1PE10.^,T32, 
1S20.6,I9,I6) 

5212   FORMAT  (//) 

213  FORMAT  (»  +  », T^-, 'POINT', T17, 'PSI*  ,T30, "VELOCITY' ,T^3» 
1' TEMPERATURE* ,T6l , 'DENSITY' ,T75, 'VISCOSITY' ,T92,'PHI' , 
2T103,'Y/A*,T116,'FQINT'  ) 

215    FORMAT  ('0*.T4,I5j1P7E15.^2X,I5) 
5211   FORMAT  ('!'  ) 

310  FORMAT  (•0»,T^,'X/A'>l?S12.6tT112,I5) 

311  FORMAT  (»+•  ,T*S«?OINT*  tT17t*  SIGMA'  ,T28  ,*  SNTHALPY*  ,  T^, 
1 • CPBAR* . T59 , * WTMIX* , T73 , ' ALPHA ( 1 ) * , T88 , *  ALPHA ( 2 ) • , 
2T103S,AL?HA(3)' ) 

312  FORMAT  (•+■ ,T103,'ALPKA(^)* ) 

a  :   ('0»,T4,I5,F13.^,1P6E15.^-) 

352  FORMAT  (I1,2I3,T10,7?9.D 

351  FORMAT  (Il,2I3,T10,^F10.Jf) 

1 971  FORMAT  ( *  +  ■ , T3 5  f  * PRAN DTL  I » ) 

3 973  FORMAT  ( •  +•  , T35 , *  MODIFIED  FERRI ' ) 

1974  FORMAT  (•+• ,T85,* MODIFIED  SCHETZ* ) 

1975  FORMAT  (•+• ,T85,' MODIFIED  ALPINIERI* ) 
350  RETURN 

END 
C 

c 


SU3R0UTINE   SOLVE 
C 
C  THIS   SUBROUTINE   SOLVES    THE  FINITE  DIFFRENCE  EQUATIONS. 

VISION    TITLE(l8),A(^),B(^),C(^),T(200),U(200), 
1WTM0L (200 ) , WTMIX ( 200 ) , RHO (200 ) , HTOT ( 200 ) , H3TAT ( 200 ) , 
2USQ (200 ) ,PSI (200 ), Y ( 200 ), EPSI ( 200 ),SG (200), CPBAR (200), 
3WU ( 200 ) ,WALFA(^,200 ) , WH(200 ) , ALFA (i+, 200 ) , HSPEC 
Wi  200),  SIGMA  (200) 
C 

COMMON    TITLE,MPSI,INDUMP,N0SPSC,NPUNCK,A,3,C,XPRTNT,XINIT 
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COMMON    X^DIAL, PRNDTL, XLEWIS,D£LPSI,XMPS,P, CI, C2fC3,C^,T,U 

C  OMMON    \L?k .WTMOL ,WTMIX , RHO , HTOT , HSTAT , HSPEG , USQ ,  PS  I , ISTS? 

C  OMMON    PCNT.DX ,  JEXIT ,Y ,EPSI , SG , CPBAR , SIGMA  ,V/U  ,WALFA  ,WH  ,NPS: 

C  OMMON    X ,  DON  IT , MHALF , IPAGS 

COMMON    PHI (200 ),3Z,C5,C6, MODEL, HEIGHT 

COMMON   XBAR(9)»IXBAR 
C 

G  DETERMINE  STEP  SIZE 

C 

C1  =1 

G2=i  ./PRNDTL 

G3=XLE,.VIS/PRNDTL 

G4= ( PRNDTL-1 . ) /PRNDTL 

G5= (XLEWIS-1 . )/PRNDTL 

FACT0R=AMAX1 ( Gl , G2 , G3 , G^ , G5 ) 

XD=FACT0R*6.*RH0(l)*U(l)*EPSI(i) 

DO   1    I=2,NPSI 

DIV=SG ( 1+1 )+SG ( I )+SG ( I )+SG( 1-1 ) 

DDX=FACT0R*DIV*1.5 

XD=AMAX1(XD,DDX) 
1  CONTINUE 

DX=  (DELPSI*DELPSI  )/(XD*XMPS ) 

DX=AMINl(DX,i|..5E-^) 

XMU=DX/DELPSI/dELPSI/XMPS 
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SOLVE  EQUATIONS  AT  THE  WALL 


WU(t)=U(l)+2.*XMU*SG(l)*(U(2)-U(l)) 

SUM=0.0 

DO  2   J=1,N0SPSC  ,  , 

WALFA(J,l)=ALFA(J,l)+XMU*SG(l)*G3*(ALFA(J,2)-ALFA(J,l))*i 

SUM=SUM+HSPEC ( J, 1 )* (ALFA (J , 2 ) -ALFA (J , 1 ) ) 
2  CGNT~ 

T1=SUM*G5 

T2=G^*(USQ(2)-USQ(1)) 

T3=G2*(KT0rT1(2)-HT0T(l) ) 

WH  ( 1  )-K  I OT  ( 1  )+2 .  *XMU  *SG  ( 1  )*  (T1+T2+T3 ) 

c 

C      SOLVE  EXPLICIT  EQUATIONS 


DO  50  I=2,NPSI 
T4=0.5*(SG(I)+SG(I+1)) 
T5=0.5>(SG(I)+SG(I-1)) 

WU(I)=U(I)+XMU*(T^*U(I+1)-(T4+T5)*U(I)+T5*U(I-1)) 
IF  (PRNDTL. GT.0.95.AND. PRNDTL. LT. 1.05)  GO  TO  20 
T6=T^*USQ(  1+1 )-  (T4+T5  )"-USQ  ( I  )+T5*USQ  ( 1-1 ) 
T66=G4*XMU*T6 
GO  TO  21 
20     T66=0.0 

IF  (XLEWIS.GT.O. 95.ANDXLEWIS.LT. 1.05)  GOTO  30 
T20=0.0 
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T21=0.0 

T22=0.0 

T23=0.0 

DO  22   J=1,NCSPSC 

T20=T20+.5*(HSPSC (J, I )+HSPEC ( J,I+1 ) )*ALFA( J.I+1 ) 

T21=T21+.5*(HSPEC(J,l)+HSPEC(J,I-l))*ALFA(J,I-l) 

T22=T22+.5*(HSPEC(J,I)+HSPEC(J,I+l))*ALFA(J,I) 

T23=T23+.5*(HSPEC(J,I)+HSPEC(J,I-1))*ALFA(JII) 
22  CONTINUE 

T202=G5*XMU*T^*T20 

T212=G5*XMU*T5*T21 

T222=G5*XMU*T^*T22 

T232=G5*XMU*T5*T23 

T77=T202+T212+T222-2.-T222-T232 

GO  TO  32 
30  T?7=0.0 

32  WH( I )=HTOT ( I )+XMU*G2* ( T4*HT0T ( 1+1 )- (T4~rT5 ) *HTOT(I )+T5* 
1HT0T(I-1))+T66+T77 

DO   33   J=1,N0SPEC 

WALFA ( J , I }=ALFA ( J , I )+G3*XMU* ( T*f*ALFA ( J , 1+1 ) - ( T4+T5 ) * 
1ALFA (J , I )+T5*ALFA ( J , 1-1 ) ) 

33  C  ON  T  PI  US 

50  CONTINUE 

DO  51    I=1,NPSI 

U(I)=WU(I) 

PHI(I)=(U(I)-U(l))/(U(MPSI)-U(l)) 

HT0T(I)=WH(I) 

USQfl)=U.(l)-;;-u(l)/50073. 

DO   51    J=1,N0SPEC 
ALFA(J,I)=WALFA(J,I) 

51  CONTINUE 
C 

G  COMPUTE   TEMPERATURE,   DENSITY,    MOLECULAR  WEIGHT,    ENTHALPIES 

G 

CALL  SUPP 
C 

DO  60   I=1,MPSI 

IF    (T(I).LS.O.O)    GO   TO   66 

DO  60   J=1,N0SPEC 

IF    (ALFA(J,I).LT.0.0.0R.ALFA(J,I).GT.1.0)    GO   TO   66 

60  CONTINUE 
GO  TO  6X 

66  X=X+DX 

WRITE    (6,67)    X.DX.ISTEP 

JEXIT=2 

PCNT=XPRINT 

GO  TO   2000 

61  IF    (JEXIT)3,^,3 
G 

C       SET  HALVING  CODE  FIRST  TIME  THRU 
C 
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>4-  JEXIT=1 

MINIT=MPSI 
MHALF=MPSI+MPSI-1 
3     X=X+DX 

IF  (MPSI-MHALF)  999,1500,1500 
G 

C      TEST  FOR  ZERO  SLOPE  AT  EDGE 
C 

999   IF  (A3S(U(N?SI)-U(I/1?SI))/U(MPSI)-0.001)  1011,1011,1004 
1011   IF  (A3S(HT0T(NPSI)-HT0T(lvIPSI))/HT0T(MPSI)-0.00l) 
1   1002,1002,1004 

1002  DO  1003  J=1,N0S?EC 

IF  ( A.3S  (ALFA  ( J ,NPSI ) -ALFA  ( J , MPSI ) ) -0 .  0005  )  100.3 ,1003 ,1004 

1003  CONTINUE 
GO  TO  2000 

C 

C      EXPAND  MESH 

C 

1004  MPSI=MPSI+1 
NPSI=MPSI-i 
GO  TO  2000 

G 

C      HALVE  MESH 
C 
1500   JSXIT=0 

DEL?SI=DSLPSI+DELPSI 

DO  ioOO  1=1, ION IT 

K=I+I-1 

U(I)=U(K) 

HT0T(I)=HT0T(K) 

USQ(I)=USQ(K) 

PHI(I)=PHI(K) 

WTMIX(I)=WTMIX(K) 

RH0(I)=RH0(K) 

T(I)=T(K) 

HSTAT ( I )=HS TAT (K) 

CP3AR(I)=CPBAR(K) 

DO  16 SO  J=1,N0SPEG 

HSPEC ( J , I ) =HS?SC ( J,K) 

ALFA(J,I)=ALFA(J,K) 
1550  CONTINUE 

XI=I-1 

PSI(I)=XI*DELPSI 
1600      CONTINUE 

MPSI=MINIT 

NPSI=MPSI-1 

DO  1700    I=MINIT,200 

U(I)  =  U(:.I?SI) 

PHI(I)=PHI(MPSI) 
HTOT(I)=HTOT(MPSI) 
V  CX    [  MIX(MPSI) 
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USQ(l)=USQ(MPSI) 

RHO(l)=RHO(MPSI) 

T(l)=T(MPSl) 

HSTAT ( I )=HSTAT(MPSI) 

CP3AR  ( I  )=CPBAR (M?SI ) 

XI=I-1 

PSI(l)=XI*DELPSI 

DO   1700   J=l,NOSPSC 

ALFA ( J , I ) =ALFA ( J ,  MPS  I ) 

HSPEG ( J ,  I )=HSPSC ( J , MPS I ) 
1700      CONTINUE 
67  FORMAT    (M"  ,T6,»NEGATIVS  SPECIES   OR  TEMPERATURE   OBTAINED.* 

!,♦    PROGRAM  TERMINATED.1 ,1P2S20.6,I10 ) 
2000      RETURN 

END 
C 
C 

C  ;- ?-'- :-— -- 

c 
c 

SUBROUTINE  SUP? 
C 

C  THIS    IS  A   COMPLEMENTARY   SUBROUTINE  TO  SUBROUTINE  SOLVE. 

G  SUP?  TAKES   VALUES    OF   TOTAL  ENTHALPY,    VELOCITY  AND   MASS 

C  FRACTIONS    FROM  SOLVE  AND   COMPUTES   STATIC    ENTHALPY,    SPECIES 

ENTHALPY,   AVERAGE  CP,    MIXTURE  MOLECULAR  WEIGHT,    DENSITY  AND 
C  TEMPERATURE  ©   EACH  PS I   GRID  POINT. 

C 

DIMENSION    TITLE(l3)tA(^),B(^),0(^),T(200),U(200), 
1WTM0L ( *0  ,WTMIX ( 200 ), RHO ( 200 ),HTOT( 200), HSTAT (200), 
2USQ(200),PSI(200),Y(200),EPSI(200),SC-(200),CP3AR(200)f 
3WU(200 )  ,WALFA (if, 200 )  ,WH(200 )  ,ALFA(4,200  ) , 
4HSPEC  (if,  200), SIGMA  (200) 
C 

COMMON   TITLE, MPSI,INDUMP,NOSPEC ,NPUNCH,A,B,C,XPRINT,XINIT 

COMMON   XFINAL,PPNDTL,XLEWIS,DELPSI,XMiPS,P,Cl,C2,C3,C^T,U 

C OMMON   ALFA ,WTMOL ,WTMIX , RHO , HTOT , HSTAT , HSPEC , U3Q , PS I , ISTE? 

COMMON    PCNT ,DX, JEXIT , Y ,EPSI ,SG,CP3AR, SIGMA ,WU ,WALFA ,WH ,NPSI 

C  OMMON   X ,  HIM  IT ,  MKALF ,  IPAGE 

COMMON    PHI ( 200 ),B 2, 05, 06 .MODEL, HEIGHT 

C OMMON   X3AR ( 9 ) , IX3AR 
C 

C  CALCULATE  MIXTURE  STATIC   ENTHALPY   AND   TEMPERATURE 

C 

DO   1    I- 1, MPS I 

HSTAT ( I )=HTOT ( I ) -USQ ( I ) 

E1=0.0 

E2=0.0 

E3=0.0 

DO   2    J=1,N0SPEC 

E1=E1+A(J)*ALFA(J,I) 
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E2=E2+B(J)*ALFA(J,I) 
E3=E3+C(J)*ALFA(J,l) 

2  CONTINUE 

IF    (E3.EQ.0.0)    GO  TO  3 

T(I)-(-E2+SQRT(E2*E2-^.*E3*(El-HSTAQ?(I))))/(2.*E3) 
GO   TO   1 

3  T(I)=(HSTAT(I)-Ei)/E2 
l          continue 

c 

C  COMPUTE  MIXTURE  MOLECULAR  WEIGHT,    CP3AR, DENSITY 

C  AND   SPECIES   ENTHALPIES 

C 

DO  5    I=1,MPSI 

sum=o.o 

DSUM=0 . 0 

DO  6   J=1,N0SPSC 

H3PEC(J,I)=A(J)+(B(J)*T(I))*T(I) 

SUM=SUM+ALFA  ( J ,  I )  /-VTMOL  ( J ) 

IF    (G(J))    40,^1,4-0 
kl  DSU  M=DSUM+B  ( J )  *ALFA  ( J ,  I ) 

GO   TO   6 
kO  D3UM=D3UM+(B(j)+2.:*C(J)*T(l)  )*ALFA(J,I) 

6  cc;jTr::;E 

CPBAR(I)=DSUM 

WTMIX(I)=SUM 

RKO ( I )= ( P*WTMIX  ( I ) )  /  ( 4-97 1 9 . 1  *T  ( I ) ) 
5  'INUE 

RETURN 

END 
C 
r 
c  __ , 

Q 

C 

SUBROUTINE  SHEAR 

C 

c  this  subroutine  models  the  turbulent  transport  mechanism 

C  PUTES   EDDY   VISCOSITY   THROUGHOUT   MPSI   GRID   POINTS. 

C 

DIMENSION    TITLE(i3),A(^),3(U),C(if),T(200),U(200). 
1WTM0L (b) ,WTMIX ( 200 ) , RHO ( 200 ) ,HTOT ( 200 ) ,H3TAT ( 200 ) , 
2U3Q(200),PSI(200),Y(200),EPSI(200),SG(200),CP3AR(200), 
3WU  ( 200  )  ,WALFA  ( k, 200  ) , WH ( 200 )  , ALFA  (4,  200  ) , 
4H3PSC (^,200), SIGMA (200) 
C 

COMMON  HITLE,MPSI,INDUMP,NOSPEC,NPUNCH,AtB,CtXPRINT,XINIT 
COMMON  XFLNAL,PP^DTL,XLET.'/IS,DEL?3I,H..:PS,P,Cl,C2,C3,Cif,T,U 
C OMMCN   ALFA , WTMGL  ,7/TMIX ,  RHO , HTOT , HSTAT  ,HSPEC  ,  U3Q , P3I ,  ISTEP 
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COMMON   PCNT,DX,JEXIT,Y,EPSI,SG,CP3AR, SIGMA, WU,WALFA,WH,NPSI 

C  OMMON   X ,  I.IIN  IT ,  MHALF ,  I  PA  GE 

COMMON    PHI (200 ) ,3Z,C5,C6, MODEL, HEIGHT 

C  OMMCN    XBAR ( 9 ) • IXBAR 
C 

C  COMPUTE  MIX  BIG   ZONE   THICKNESS 

C 

IF    (U(MPSl).LT.U(l))    GO  TO  5 

WUMAX=.95*U(MPSI) 

WUMLN=1.05U(1) 

DO   1    I=1,MPSI 

IMAX=I 

IF    (U(I).GE.WUMAX)    GO   TO   2 

1  CONTINUE 

2  CONTINUE 

DO   3    I=1,MPSI 

IMIN=I 

IF    (U(I).GE.WUMIN)    GO   TO    >4 

3  CONTINUE 

5  WUMAX=1.^05*U(MPSI) 
WUMIN=.95*U(1) 

DO   6    1=1, MPS I 

IMAX=I 

IF    (U(I).LE.WUMAX)    GO  TO  7 

6  CONTINUE 

7  CONTINUE 

DO  8    1=1, MPS I 

IMIN=I 

IF    (U(I).LE.WUMIN)    GO   TO   k 

8  CONTINUE 
k  CONTINUE 

DUM2=WUMIN-U(IMIN) 
DUM3=U ( IMIN-1 )  -U  ( IMIN  ) 
DUM4=Y  ( IMIN-1  )-Y  ( IMIN  ) 
YMIN=DUM2*DUMVDUM3+Y  ( IMIN  ) 
DUM5=WUT/tAX-U  ( IMAX ) 

dum6=u  ( imax-i )  -u  (  imax  ) 

DUM7=Y  ( IMAX-1 )  -Y  ( IMAX ) 

YMAX=DUM5*DUM7/DUM6+Y  ( IMAX ) 

BZ=YMAX-YMIN 

RETURN 
C 

ENTRY   SHEAR1    (*) 
C 

C  MODEL   USED    :    ************   PRANDTL    I    ************ 

C 

EPSI(1)=0.0 

DO   11    I=2,MPSI 

DFDY=A33 ( ( U ( 1+1 )-U( 1-1 ) )/(2 . *DELPSI )  ) 

EPSI ( I )=C1*C1*BZ*BZ*RH0( I )*RHO( I )*U ( l)*DYDY 
11  CONTINUE 
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99  RETURN  1 

100  RETURN 
C 

RY  SHEAR3  (*) 
C 

C      MODEL  USED  »  ************  MODIFIED  FERRI  ************ 
C 

RUMAX=RHO(l)*U(l) 

RUMIN=RHO(i)*U(l) 

DO  13  I=2,MPSI 

RUX=RHO(l)*U(I) 

RUMAX=AMAXI  ( RUMAX ,  RUX ) 

RUMLN=AMIN 1 ( RUMLN , RUX ) 

13  CONTINUE 
DELRU=RUMAX-RUMIN 
VISC=C3*BZ*DELRU 
DO   Ik   1=1, MPS I 
EPSI(I)=VISC 

14  CONTINUE 
11Q  RETURN  1 
120        RETURN 

C 

ENTRY  SHEARS  (*) 
C 

G      MODEL  USED  i    ************  SCHETZ  ************ 
G 

DUM=0.0 

DO  15  I=1,MPSI 

DUM=DUM+ABS((RH0(MPSI)*U(MPSI)/RH0(I)/U(I))-1.) 

15  CONTINUE 
VISC=DUM*DELPS  I  *C  k 
DO   21   I=1,MPSI 
EPSI(I)=VISC 

21  CONTINUE 

129  RETURN    1 

130  RETURN 
C 

ENTRY  SHEAR5    (*) 
C 

C  MODEL  USED    »    ************  MODIFIED  ALPINIERI   ************ 

C 

RATRO=RHO ( MPS I )/RHO ( 1 ) 

RATMOM=(RHO(MPSI)nj(MPSl)*U(MPSl))/(RHO(l)*U(i)*U(i)) 

VISC=C5->32*RH0(1)^U(1)*(RATR0+RA^: 

DO   16    1=1,  LIPS  I 

EPSI(I)=VISC 

16  CONTINUE 
139  RETURN  1 
1^0        RE 1 

END 
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